Controller for controlling a plant

ABSTRACT

The present invention provides a controller for controlling a modeled plant robustly against disturbance. The controller comprises an estimator and a control unit. The estimator estimates disturbance applied to the plant. The control unit determines an input to the plant so that an output of the plant converges to a desired value. The input to the plant is determined to include a value obtained by multiplying the estimated disturbance by a predetermined gain. Since estimated disturbance is reflected in the input to the plant, control having robustness against disturbance is implemented. The controller may comprise a state predictor. The state predictor predicts the output of the plant based on the estimated disturbance and dead time included in the plant. The control unit determined the input to the plant so that the predicted output converges to a desired value. Since the state predictor allows for the dead time, the accuracy of the control is improved. The estimated disturbance is reflected in the predicted output, an error between the predicted output and an actual output of the plant is removed.

BACKGROUND OF THE INVENTION

[0001] The present invention relates to an apparatus for robustlycontrolling a plant against disturbance.

[0002] The amount of air introduced into an engine is typicallycontrolled so as to achieve a desired engine torque. According to aconventional method, a desired amount of air introduced into the engineis determined referring to a map based on an opening angle of anaccelerator pedal, a vehicle speed and a selected transmission gearratio. An opening angle of a throttle valve is controlled in accordancewith the desired amount of air to be introduced into the engine.

[0003] According to another method disclosed in Japanese Patent No.2780345, a desired torque of an engine output shaft is determined inaccordance with an opening angle of an accelerator pedal and arotational speed of the output shaft of a torque converter. A desiredopening angle of a throttle valve is determined referring to apredetermined table based on an actual rotational speed and the desiredtorque of the engine output shaft. Air is introduced into the engine inaccordance with the desired throttle opening angle.

[0004] Such conventional control does not take into account disturbanceapplied into the intake manifold and dead time from the throttle valveto the engine. Such factors reduce the accuracy of controlling air to beintroduced into the engine, causing vibration in the engine torque.

[0005] Such a problem regarding robustness against disturbance exists inrotational speed control for the engine.

[0006] Conventionally, when an engine is idling, a conventional PIDcontrol is performed for controlling an engine rotational speed.According to such conventional rotational speed control, when a suddenchange in the engine load occurs during idling operation, the enginetends to stop because the engine rotational speed cannot be stable. Forexample, when a vehicle with a manual transmission mechanism starts anda clutch is forced to engage suddenly, the engine of the vehicle tendsto stop.

[0007] Japanese Patent No. 3203602 discloses a scheme for reducing shockthat may make passengers of the vehicle uncomfortable when gear changeoccurs in an automatic transmission mechanism. According to the scheme,a desired torque of a driving shaft is determined based on a vehiclespeed and an opening angle of the accelerator pedal. A desired enginetorque and a desired engine rotational speed that provide the desireddriving shaft torque are determined. An opening angle of the throttlevalve is determined based on the desired engine torque and the desiredengine rotational speed. The throttle valve is controlled in accordancewith the determined opening angle.

[0008] Conventionally, in a vehicle that comprises an automatic manualtransmission (automatic MT) or an automatic transmission (AT), when gearchange occurs, a rotational speed synchronization control capable ofachieving a quick response is not performed. Thus, a rotational speedcannot be adapted to a selected transmission gear ratio quickly.

[0009] Thus, there is a need for control that has high robustnessagainst disturbance.

SUMMARY OF THE INVENTION

[0010] According to one aspect of the present invention, a controllerfor controlling a modeled plant is provided. The controller comprises anestimator for estimating disturbance applied to the plant and a controlunit for determining an input to the plant so that an output from theplant converges to a desired value. The input to the plant is determinedto include a value obtained by multiplying the estimated disturbance bya predetermined gain.

[0011] An error between the output of the plant and a desired value forthe output of the plant may be caused by disturbance applied to theplant. According to the invention, the controller can quickly cause suchan error to converge to zero since the input to the plant includes theestimated disturbance.

[0012] According to one embodiment of the present invention, the controlunit uses a preview control to determine the input to the plant. Thepreview control can implement a feasible control for the plant when deadtime is included in the plant.

[0013] According to one embodiment of the present invention, the inputto the plant is determined to include a value obtained by multiplying adesired value for the output of the plant by a predetermined gain. Thus,the capability that the output of the plant follows the desired value isimproved.

[0014] According to another embodiment of the present invention, thecontrol unit uses a response assignment control to determine the inputto the plant. The response assignment control enables an error betweenthe output of the plant and a desired value for the output of the plantto converge to zero without generating overshooting.

[0015] According to one embodiment of the present invention, theestimator is an adaptive disturbance observer that identifies thedisturbance using a recursive identification algorithm. Such a recursiveidentification algorithm can quickly and stably identify the estimateddisturbance. When noise is included in the output of the plant,variation may occur in the estimated disturbance due to such noise. Theeffect of a statistical process of the recursive identificationalgorithm can remove such variation in the estimated disturbance.

[0016] According to another aspect of the present invention, thecontroller further comprises a state predictor for predicting the outputof the plant based on the estimated disturbance and dead time includedin the plant. The control unit determines the input to the plant so thatthe predicted output converges to a desired value for the output of theplant. Since the dead time is compensated for by the state predictor, aresponse of the control is improved. Since the predicted output isdetermined taking into account the estimated disturbance, an errorbetween the predicted output and the actual output of the plant isremoved.

[0017] Conventional generalized predictive control requires decreasing again when a dead time of the plant is taken into account. The statepredictor removes such decrease of the gain since the dead time iscompensated for by the state predictor.

[0018] According to one embodiment of the present invention, the plantis an intake manifold connected to an engine. The intake manifold ismodeled so that its input is a desired value for an opening angle of avalve that controls an amount of air introduced into the intake manifoldand its output is an amount of air introduced into the engine. Thus, anamount of air introduced into the engine converges to a desired valuewith high accuracy, thereby accurately controlling an engine torque. Theinput into the plant may be a desired value for an opening angle of athrottle valve provided in the intake manifold.

[0019] According to one embodiment of the present invention, a modelparameter for the modeled plant is determined based on an actual enginerotational speed and an actual opening angle of the throttle valve. Themodel parameter thus determined achieves an accurate control for theengine torque under various engine operating conditions.

[0020] According to one embodiment of the present invention, the plantis an engine. The engine is modeled so that its input is a desired valuefor an amount of air introduced into the engine and its output is arotational speed of the engine. Thus, engine stall that may occur whenthe engine starts is suppressed. A response of the engine rotationalspeed control when a transmission gear change occurs is improved.

[0021] According to another embodiment of the present invention, thecontroller determines a model parameter for the modeled plant based on adetected rotational speed. The input to the plant is determined usingthe model parameter. The model parameter thus determined achieves anaccurate control for the rotational speed under various engine operatingconditions.

[0022] According to another embodiment of the present invention, theinput to the plant includes a value obtained by multiplying by apredetermined gain an estimated value for a torque required for drivingthe vehicle. According to another embodiment of the present invention,the input to the plant includes a value obtained by multiplying by apredetermined gain an estimated value for a torque required for drivingequipments mounted on the vehicle. Thus, an error between the output ofthe plant and its desired value that may be caused by thevehicle-driving torque and the equipment-driving torque can converge.

[0023] According to one embodiment of the present invention, the statepredictor further determines the predicted output based on the estimatedvalue for the vehicle-driving torque. According to another embodiment ofthe present invention, the state predictor determines the predictedoutput based on the estimated value for the equipment-driving torque.Thus, an error between the predicted output and a desired value for theoutput of the plant that may be caused by the vehicle-driving torque andthe equipment-driving torque can converge.

BRIEF DESCRIPTION OF THE DRAWINGS

[0024]FIG. 1 is a block diagram of an internal combustion engine and itscontrol unit in accordance with one embodiment of the present invention.

[0025]FIG. 2 is a block diagram of a control unit in accordance with oneembodiment of the present invention.

[0026]FIG. 3 shows a structure of an intake air amount feedback controlin accordance with one embodiment of the present invention.

[0027]FIG. 4 shows a virtual controlled object in a simulation for anintake air amount feedback control in accordance with one embodiment ofthe present invention.

[0028]FIG. 5 shows a result of a case G-1 for an intake air amountcontrol simulation in accordance with one embodiment of the presentinvention.

[0029]FIG. 6 shows a result of a case G-2 for an intake air amountcontrol simulation in accordance with one embodiment of the presentinvention.

[0030]FIG. 7 shows a result of a case G-3 for an intake air amountcontrol simulation in accordance with one embodiment of the presentinvention.

[0031]FIG. 8 shows a result of a case G-4for an intake air amountcontrol simulation in accordance with one embodiment of the presentinvention.

[0032]FIG. 9 shows a result of a case G-5 for an intake air amountcontrol simulation in accordance with one embodiment of the presentinvention.

[0033]FIG. 10 shows a result of a case G-6 for an intake air amountcontrol simulation in accordance with one embodiment of the presentinvention.

[0034]FIG. 11 shows a result of a case G-7 for an intake air amountcontrol simulation in accordance with one embodiment of the presentinvention.

[0035]FIG. 12 shows a comparison between the case G-8 and the case G-7for the intake air amount control simulation in accordance with oneembodiment of the present invention.

[0036]FIG. 13 shows a structure of a rotational speed feedback controlin accordance with one embodiment of the present invention.

[0037]FIG. 14 shows a switching line of a response assignment control inaccordance with one embodiment of the present invention.

[0038]FIG. 15 shows a relation between a convergence speed and the valueof a setting parameter of a switching function for a response assignmentcontrol in accordance with one embodiment of the present invention.

[0039]FIG. 16 shows a virtual controlled object in a simulation for arotational speed feedback control in accordance with one embodiment ofthe present invention.

[0040]FIG. 17 shows a result of a case N-1 for a rotational speedcontrol simulation in accordance with one embodiment of the presentinvention.

[0041]FIG. 18 shows a result of a case N-2 for a rotational speedcontrol simulation in accordance with one embodiment of the presentinvention.

[0042]FIG. 19 shows a result of a case N-3 for a rotational speedcontrol simulation in accordance with one embodiment of the presentinvention.

[0043]FIG. 20 shows a result of a case N-4 for a rotational speedcontrol simulation in accordance with one embodiment of the presentinvention.

[0044]FIG. 21 shows a result of a case N-5 for a rotational speedcontrol simulation in accordance with one embodiment of the presentinvention.

[0045]FIG. 22 shows a result of a case N-6 for a rotational speedcontrol simulation in accordance with one embodiment of the presentinvention.

[0046]FIG. 23 shows a result of a case N-7 for a rotational speedcontrol simulation in accordance with one embodiment of the presentinvention.

[0047]FIG. 24 shows a result of a case N-8 for a rotational speedcontrol simulation in accordance with one embodiment of the presentinvention.

[0048]FIG. 25 shows a result of a case N-9 for a rotational speedcontrol simulation in accordance with one embodiment of the presentinvention.

[0049]FIG. 26 shows a comparison between the case N-8 and the case N-9for the rotational speed control simulation in accordance with oneembodiment of the present invention.

[0050]FIG. 27 shows a flowchart of a rotational speed feedback controland an intake air amount feedback control in accordance with oneembodiment of the present invention.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0051] Structure of Internal Combustion Engine and Control Unit

[0052] Referring to the drawings, specific embodiments of the inventionwill be described. FIG. 1 is a block diagram showing an internalcombustion engine (hereinafter referred to as an engine) and its controlunit in accordance with one embodiment of the invention.

[0053] An electronic control unit (hereinafter referred to as an ECU) 1comprises an input interface 1 a for receiving data sent from each partof the vehicle, a CPU 1 b for carrying out operations for controllingeach part of the vehicle, a memory 1 c including a read only memory(ROM) and a random access memory (RAM), and an output interface 1 d forsending control signals to each part of the vehicle. Programs andvarious data for controlling each part of the vehicle are stored in theROM. The ROM may be a rewritable ROM such as an EEPROM. The RAM provideswork areas for operations by the CPU 1 b, in which data sent from eachpart of the vehicle as well as control signals to be sent out to eachpart of the vehicle are temporarily stored.

[0054] The engine 2 is, for example, an engine equipped with fourcylinders. Each cylinder comprises an intake valve 5 for connecting acombustion chamber 7 to an intake manifold 3 and an exhaust valve 6 forconnecting the combustion chamber 7 to an exhaust manifold 4.

[0055] A throttle valve 8 is disposed upstream of the intake manifold 3.The throttle valve 8 is an electronic control valve. An opening angle ofthe throttle valve 8 is controlled by a control signal from the ECU 1. Athrottle valve opening (θ TH) sensor 9, which is connected to thethrottle valve 8, outputs an electric signal corresponding to an openingangle of the throttle valve 8 and sends the electric signal to the ECU1.

[0056] An airflow meter (AFM) 10 is provided upstream of the throttlevalve 8. The airflow meter 10 detects the amount Gth of air passingthrough the throttle valve 8, and sends it to the ECU 1. The airflowmeter 10 may be a vane-type airflow meter, a Karman-vortex type airflowmeter, a hot-wire type airflow meter or the like.

[0057] An intake manifold pressure (Pb) sensor 11 is provided in theintake manifold 3 downstream of the throttle valve 8. A pressure Pb ofthe intake manifold detected by the Pb sensor 11 is sent to the ECU 1.

[0058] A fuel injection valve 12 is installed for each cylinder in theintake manifold 3 upstream of the intake valve 5. The fuel injectionvalve 12 is supplied with fuel from a fuel tank (not shown) and drivenin accordance with a control signal from the ECU 1.

[0059] A rotational speed (Ne) sensor 13 is attached to the periphery ofthe camshaft or the periphery of the crankshaft (not shown) of theengine 2, and outputs a CRK signal at a predetermined crank angle cycle(for example, a cycle of 30 degrees). The cycle length of the CRK signalis shorter than the cycle length of a TDC signal that is issued at acrank angle cycle associated with a TDC position of the piston. Pulsesof the CRK signal are counted by the ECU 1 to determine the rotationalspeed Ne of the engine.

[0060] Signals sent to the ECU 1 are passed to the input interface 1 a.The input interface 1 a converts analog signal values into digitalsignal values. The CPU 1 b processes the resulting digital signals,performs operations in accordance with the programs stored in the ROM 1c, and creates control signals. The output interface 1 d sends thesecontrol signals to actuators for the fuel injection valve 12 and otheractuators.

[0061] Air introduced into the intake manifold 3 through the throttlevalve 8 is filled in the chamber 14. When the intake valve 5 is opened,the air in the chamber 14 is supplied to the combustion chamber 7 of theengine 2. Fuel is supplied from the fuel injection valve 12 to thecombustion chamber 7. The air-fuel mixture is ignited by a spark plug(not shown) in the combustion chamber 7.

[0062] Block Diagram of Control Unit

[0063]FIG. 2 shows a block diagram of the control unit in accordancewith one embodiment of the present invention. An engine torque settingunit 20, a rotational speed feedback (FB) controller 21, a switch 22 andan intake air amount feedback (FB) controller 23 are typicallyimplemented by computer programs that are stored in the memory 1 c (FIG.1). Alternatively, functions of the blocks may be implemented bysoftware, firmware, hardware or any combination thereof.

[0064] The torque setting unit 20 refers to a map pre-stored in thememory 1 c based on an opening angle of the accelerator pedal, a vehiclespeed, a transmission gear ratio or the like to determine a desiredengine torque. The torque setting unit 20 determines the amount ofintake air required for the desired engine torque as a desired intakeair amount Gcyl_cmd.

[0065] The rotational speed FB controller 21 performs a feedback controlfor the engine rotational speed NE. The engine 2 is an object to becontrolled (hereinafter referred to as a “plant”) by the rotationalspeed feedback control. In the feedback control, a control input is adesired intake air amount Gcyl_cmd and a control output is the enginerotational speed NE. The rotational speed FB controller 21 determinesthe desired intake air amount Gcyl_cmd so that the engine rotationalspeed NE converges to a desired value.

[0066] When the vehicle is in a normal running condition, an intake airamount FB controller 23 is connected to the torque setting unit 20through the switch 22. The intake air amount FB controller 23 uses thedesired intake air amount Gcyl_cmd determined by the torque setting unit20. When the engine is idling or when gear change is being carried outin the transmission, the intake air amount FB controller 23 is connectedto the rotational speed FB controller 21. The intake air amount FBcontroller 23 uses the desired intake air amount Gcyl_cmd determined bythe rotational speed FB controller 21.

[0067] The intake air amount FB controller 23 performs a feedbackcontrol for the amount of intake air Gcyl introduced into the cylinderof the engine. The plant of the feedback control is the intake manifold3. In the feedback control, a control input is a desired value THcmd forthe opening angle of the throttle valve and a control output is theintake air amount Gcyl introduced into the engine. The intake air amountFB controller 23 determines the desired value THcmd for the openingangle of the throttle valve so that the intake air amount Gcyleconverges to the desired value Gcyl_cmd. The throttle valve 8 iscontrolled by the ECU 1 in accordance with the desired throttle openingangle THcmd.

[0068] Thus, when the engine is idling or when gear change occurs in thetransmission, the intake air amount for causing the rotational speed NEto converge to a desired value is established as a desired intake airamount. Accordingly, engine stall when the engine is idling can besuppressed. The rotational speed when gear change occurs in thetransmission can stably and quickly converge to a desired value.

[0069] In the present specification, the intake air amount FB controlwill be first described and then the rotational speed FB control will bedescribed.

[0070] 1. Intake Air Amount Feedback Control

[0071] 1.1 Modeling of Dynamic Behavior of Intake Air

[0072] A method for modeling the dynamic behavior of intake air will bedescribed. The intake manifold 3 is represented by a model in which itsinput is THcmd and its output is Gcyl.

[0073] The amount of intake air Gcyl′ introduced into each cylinder ineach cycle can be expressed by equation (1) based on the ideal gasequation of state that is known. In the equation (1), Kηc′ denotes acharging efficiency (%) of the intake Omanifold, Pb denotes a pressure(Pa) of the intake manifold, Vcyl denotes a volume (m³) of the cylinder,Tcyl denotes a temperature inside of the cylinder (K), R denotes the gasconstant (m³▪Pa/g▪K), and “n” denotes an identifier for identifying eachsampling cycle. $\begin{matrix}{{{{Gcy1}^{\prime}(n)} = {K\quad \eta \quad {c \cdot {{Pb}(n)}}}}{{{where}\quad K\quad \eta \quad c} = {K\quad \eta^{\prime}c\frac{Vcy1}{{R \cdot T}\quad {cy1}}}}} & (1)\end{matrix}$

[0074] In case of an in-line 4-cylinder engine, the intake of air isperformed twice for each rotation of the engine. The amount of airintroduced into the cylinder per unit time Gcyl is shown in the equation(2). Here, NE denotes an engine rotational speed (rpm) and k denotes anidentifier for identifying each sampling cycle. Fcyl is a function ofthe rotational speed NE. $\begin{matrix}\begin{matrix}{{{Gcy1}(k)} = \frac{{{Gcy1}^{\prime}(k)} \cdot 2 \cdot {NE}}{60}} \\{\quad {= {{Fcy1} \cdot {{Pb}(k)}}}} \\{{{where}\quad {Fcy1}} = \frac{{2 \cdot K}\quad \eta \quad {c \cdot N}\quad E}{60}}\end{matrix} & (2)\end{matrix}$

[0075] On the other hand, the amount of air ΔGb that is to be filled inthe chamber 14 is shown by the equation (3).

ΔGb(k)=Gth(k)−Gcyl(k)  (3)

[0076] As to the chamber 14, the equation (5) is derived from the idealgas equation of state (4). Pb, Vb and Tb denote a pressure (Pa), avolume (m³) and a temperature (K) of the intake manifold, respectively.R indicates the gas constant as described above.

Pb(k)·Vb=Gb(k)·R·Tb  (4)

[0077]

ΔPb(k)·Vb=Gb(k)·R·Tb

[0078]

$\begin{matrix}{{\Delta \quad G\quad {b(k)}} = \frac{\Delta \quad P\quad {{b(k)} \cdot V}\quad b}{{R \cdot T}\quad b}} & (5)\end{matrix}$

[0079] The equation (6) is obtained by substituting the equation (5)into the equation (3). The amount of intake air Gcyl is represented as afunction of the pressure Pb of the intake manifold as shown by theequation (6). T denotes the length of the sampling cycle.$\begin{matrix}{\frac{\Delta \quad P\quad {{b(k)} \cdot V}\quad b}{{R \cdot T}\quad b} = {{G\quad t\quad {h(k)}} - {G\quad c\quad {{y1}(k)}}}} & (6) \\\left. \quad\Downarrow \right. & \quad \\{\frac{\frac{\left( {{{Pb}(k)} - {{Pb}\left( {k - 1} \right)}} \right)}{T}{Vb}}{R \cdot {Tb}} = {{{Gth}(k)} - {G\quad {{cy1}(k)}}}} & \quad \\\left. \quad\Downarrow \right. & \quad \\{{G\quad {{cy1}(k)}} = {{\frac{- {Vb}}{R \cdot {Tb} \cdot T}{{Pb}(k)}} + {\frac{Vb}{R \cdot {Tb} \cdot T}{{Pb}\left( {k - 1} \right)}} + {{Gth}(k)}}} & \quad\end{matrix}$

[0080] In order to use Gcyl to express the Pb of the equation (6), theequation (7) is derived by substituting the equation (2) into theequation (6). The equation (7) represents a model of the dynamicbehavior of intake air in which its input model is Gth. $\begin{matrix}\begin{matrix}{{{G\quad {{cy1}(k)}} = {{\frac{- {Vb}}{{Fcy1} \cdot R \cdot {Tb} \cdot T}G\quad {{cy1}(k)}} +}}\quad} \\{\quad {{\frac{Vb}{{Fcy1} \cdot R \cdot {Tb} \cdot T}G\quad {{cy1}\left( {k - 1} \right)}} + {{Gth}(k)}}\quad} \\\left. \Downarrow\quad \right. \\{{{G\quad {{cy1}(k)}} = {{\frac{Vb}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}G\quad {{cy1}\left( {k - 1} \right)}} +}}\quad} \\{\quad {\frac{{Fcy1} \cdot R \cdot {Tb} \cdot T}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}{{Gth}(k)}}\quad}\end{matrix} & (7)\end{matrix}$

[0081] On the other hand, a relationship between the amount of air Gthpassing through the throttle valve and an opening angle TH of thethrottle valve is expressed by the equation (8). Here, Pc represents apressure upstream of the throttle valve. Fth denotes a flow rate pereffective opening angle of the throttle valve (g/deg), which isdetermined in accordance with the pressure Pb downstream of the throttlevalve (that is, the pressure of the intake manifold) and the pressure Pcupstream of the throttle valve. The equation (9) is obtained bysubstituting the equation (8) into the equation (7). The equation (9)represents a model of the dynamic behavior of intake air in which itsinput is the opening angle TH of the throttle valve.

Gth(k)=Fth·TH(k)  (8)

[0082] $\begin{matrix}\begin{matrix}{{{G\quad {{cy1}(k)}} = {{\frac{Vb}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}G\quad {{cy1}\left( {k - 1} \right)}} +}}\quad} \\{\quad {\frac{{Fth} \cdot {Fcy1} \cdot R \cdot {Tb} \cdot T}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}{{TH}(k)}}\quad}\end{matrix} & (9)\end{matrix}$

[0083] A relationship between a desired throttle opening angle THcmd andthe actual throttle opening angle TH of the electronic throttle valve isrepresented by the equation (10). The equation (10) is a first-orderdelay system having a dead time “dth.” The dead time dth is mainlycaused by electronic communication that is required for operating thethrottle valve. The equation (11) is obtained by substituting theequation (10) into the equation (9).

TH(k)=Ath·TH(k−1)+(1−Ath)·THcmd(k−dth)  (10)

[0084] $\begin{matrix}{{G\quad {{cy1}(k)}} = {{\frac{Vb}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}G\quad {{cy1}\left( {k - 1} \right)}} +}} & (11) \\{\quad {{\frac{{Ath} \cdot {Fth} \cdot {{Fcy1}({NE})} \cdot R \cdot {Tb} \cdot T}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}{{TH}\left( {k - 1} \right)}}\quad +}} & \quad \\{\quad {\frac{\left( {1 - {Ath}} \right) \cdot {Fth} \cdot {Fcy1} \cdot R \cdot {Tb} \cdot T}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}{{THcmd}\left( {k - {dth}} \right)}}\quad} & \quad\end{matrix}$

[0085] It is seen from the equation (9) that TH(k−1) can be expressedusing Gcyl(k−1) and Gcyl(k−2). The equation (12) is obtained bysubstituting the TH(k−1) into the equation (11). The equation (12) is amodel equation of the dynamic behavior of intake air in which its inputis the desired throttle opening angle THcmd and its output is the amountof intake air Gcyl. $\begin{matrix}{{G\quad {{cy1}(k)}} = {{\frac{Vb}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}G\quad {{cy1}\left( {k - 1} \right)}} +}} & (12) \\{\quad {{{Ath} \cdot {{Gcy1}\left( {k - 1} \right)}} - {\frac{{Ath} \cdot {Vb}}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}G\quad {{cy1}\left( {k - 2} \right)}} +}} & \quad \\{\quad {\frac{\left( {1 - {Ath}} \right) \cdot {Fth} \cdot {Fcy1} \cdot R \cdot {Tb} \cdot T}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}{{THcmd}\left( {k - {dth}} \right)}}\quad} & \quad \\\left. \quad\Downarrow \right. & \quad \\{{{Gcy1}(k)} = {{A\quad {{air1} \cdot {{Gcy1}\left( {k - 1} \right)}}} + {A\quad {{air2} \cdot {{Gcy1}\left( {k - 2} \right)}}} +}} & \quad \\{\quad {{Bair1} \cdot {{THcmd}\left( {k - {dth}} \right)}}} & \quad \\{{{where}\quad {Aair1}} = \frac{{{Ath}\left( {{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}} \right)} + {Vb}}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}} & \quad \\{\quad {{Aair2} = {- \frac{{Ath} \cdot {Vb}}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}}}} & \quad \\{\quad {{Bair1} = \frac{\left( {1 - {Ath}} \right) \cdot {Fth} \cdot {Fcy1} \cdot R \cdot {Tb} \cdot T}{{{Fcy1} \cdot R \cdot {Tb} \cdot T} + {Vb}}}} & \quad\end{matrix}$

[0086] Model parameters Aair1, Aair2 and Bair1 include Fcyl and Fth thatvary in accordance with the rotational speed NE, the intake manifoldpressure Pb and the pressure Pc upstream of the throttle valve. Themodel parameters corresponding to the rotational speed NE and thethrottle opening angle TH may be pre-stored in the memory lc as a map.Alternatively, the controller may comprise an identifier for identifyingsuch model parameters.

[0087] 1.2 Problem in Applying a Generalized Predictive Control (GPC)

[0088] According to the invention, feedback control for the amount ofintake air is implemented by a preview control algorithm. As a schemesimilar to the preview control, generalized predictive control(hereinafter referred to as GPC) is known (in some cases, the GPC isincluded in the category of preview control). However, it is impossibleto construct a feasible intake air amount feedback controller 23 only byusing such conventional GPC. The reason will be described.

[0089] The model for the dynamic behavior of intake air shown in theequation (12) can be expressed as shown by the equation (13). Here, itis assumed that a value of the dead time dth is “2”.

Gcyl(k)=Aair1·Gcyl(k−1)+Aair2·Gcyl(k−2)+Bair1·THcmd(k−2)

Gcyl(k+1)=Aair1·Gcyl(k)+Aair2·Gcyl(k−1)+Bair1·THcmd(k−1)  (13)

[0090] When the equation (13) is expressed by a state-space equation,the equation (14) is obtained. $\begin{matrix}{\begin{bmatrix}\begin{matrix}{{Gcy1}\left( {k + 1} \right)} \\{{Gcy1}(k)}\end{matrix} \\{{THcmd}(k)}\end{bmatrix} = {{\begin{bmatrix}{Aair1} & {Aair2} & {Bair1} \\1 & 0 & 0 \\0 & 0 & 0\end{bmatrix}\begin{bmatrix}\begin{matrix}{{Gcy1}(k)} \\{{Gcy1}\left( {k - 1} \right)}\end{matrix} \\{{THcmd}\left( {k - 1} \right)}\end{bmatrix}} +}} & (14) \\{\quad {\begin{bmatrix}\begin{matrix}0 \\0\end{matrix} \\1\end{bmatrix}{{THcmd}(k)}}} & \quad \\{{{Gcy1}(k)} = {\begin{bmatrix}1 & 0 & 0\end{bmatrix}\begin{bmatrix}\begin{matrix}{{Gcy1}(k)} \\{{Gcy1}\left( {k - 1} \right)}\end{matrix} \\{{THcmd}\left( {k - 1} \right)}\end{bmatrix}}} & \quad \\\left. \quad\Downarrow \right. & \quad \\{{X^{\prime}\left( {k + 1} \right)} = {{\Phi^{\prime}{X^{\prime}(k)}} + {G^{\prime}{{THcmd}(k)}}}} & \quad \\{{{Gcy1}(k)} = {C^{\prime}{X^{\prime}(k)}}} & \quad \\{{{where}\quad \Phi^{\prime}} = \begin{bmatrix}{Aair1} & {Aair2} & {Bair1} \\1 & 0 & 0 \\0 & 0 & 0\end{bmatrix}} & \quad \\{G^{\prime} = {{\begin{bmatrix}\begin{matrix}0 \\0\end{matrix} \\1\end{bmatrix}\quad C^{\prime}} = {{\begin{bmatrix}1 & 0 & 0\end{bmatrix}\quad {X^{\prime}(k)}} = \begin{bmatrix}\begin{matrix}{{Gcy1}(k)} \\{{Gcy1}\left( {k - 1} \right)}\end{matrix} \\{{THcmd}\left( {k - 1} \right)}\end{bmatrix}}}} & \quad\end{matrix}$

[0091] A difference operator Δ that is defined as Δ=1−Z⁻¹ is introducedto define an augmented system as shown by the equation (15). In theaugmented system, one additional row is introduced so as to derive anintegral term for suppressing a steady-state error. $\begin{matrix}{\begin{bmatrix}\begin{matrix}{{Gcy1}\left( {k + 1} \right)} \\{\Delta \quad {{Gcy1}(k)}} \\{\Delta \quad {{Gcy1}\left( {k - 1} \right)}}\end{matrix} \\{\Delta \quad {{THcmd}(k)}}\end{bmatrix} = {{\begin{bmatrix}1 & {Aair1} & {Aair2} & {Bair1} \\0 & {Aair1} & {Aair2} & {Bair1} \\0 & 1 & 0 & 0 \\0 & 0 & 0 & 0\end{bmatrix}\begin{bmatrix}\begin{matrix}{{Gcy1}(k)} \\{\Delta \quad {{Gcy1}\left( {k - 1} \right)}} \\{\Delta \quad {{Gcy1}\left( {k - 2} \right)}}\end{matrix} \\{\Delta \quad {{THcmd}\left( {k - 1} \right)}}\end{bmatrix}} +}} & (15) \\{\quad {\begin{bmatrix}\begin{matrix}0 \\0 \\0\end{matrix} \\1\end{bmatrix}{{THcmd}(k)}}} & \quad \\{{{Gcy1}(k)} = {\begin{bmatrix}1 & 0 & 0 & 0\end{bmatrix}\begin{bmatrix}\begin{matrix}{{Gcy1}(k)} \\{\Delta \quad {{Gcy1}\left( {k - 1} \right)}} \\{\Delta \quad {{Gcy1}\left( {k - 2} \right)}}\end{matrix} \\{\Delta \quad {{THcmd}\left( {k - 1} \right)}}\end{bmatrix}}} & \quad \\\left. \quad\Downarrow \right. & \quad \\{{X^{\prime}\left( {k + 1} \right)} = {{\Phi \quad {X^{\prime}(k)}} + {{G \cdot \Delta}\quad {{THcmd}(k)}}}} & \quad \\{{{Gcy1}(k)} = {C \cdot {X^{\prime}(k)}}} & \quad \\{\quad {\Phi = \begin{bmatrix}1 & {Aair1} & {Aair2} & {Bair1} \\0 & {Aair1} & {Aair2} & {Bair1} \\0 & 1 & 0 & 0 \\0 & 0 & 0 & 0\end{bmatrix}}} & \quad \\{G = {{\begin{bmatrix}\begin{matrix}0 \\0 \\0\end{matrix} \\1\end{bmatrix}\quad C} = {{\begin{bmatrix}1 & 0 & 0 & 0\end{bmatrix}\quad {X^{\prime}(k)}} = \begin{bmatrix}\begin{matrix}{{Gcy1}(k)} \\{\Delta \quad {{Gcy1}\left( {k - 1} \right)}} \\{\Delta \quad {{Gcy1}\left( {k - 2} \right)}}\end{matrix} \\{\Delta \quad {{THcmd}\left( {k - 1} \right)}}\end{bmatrix}}}} & \quad\end{matrix}$

[0092] The GPC is a technique of causing the controlled variable Gcyl toconverge to the desired value Gcyl_cmd in a time period M from a time(k) to a time (k+M). A cost function JG as shown by the equation (16) isdefined, where H is a weighting parameter (>0). $\begin{matrix}{J_{G} = {\sum\limits_{j = 1}^{M}\quad \left\lbrack {\left\{ {{{Gcyl\_ cmd}\left( {k + j} \right)} - {{Gcyl}\left( {k + j} \right)}} \right\}^{2} + {H\quad \Delta \quad {{THcmd}\left( {k + j - 1} \right)}^{2}}} \right\rbrack}} & (16)\end{matrix}$

[0093] A control input ΔTHcmd that minimize the cost function J_(G) canbe determined by using the principle of optimality. The control input ΔTHcmd is expressed as shown by the equation (17) by using solution P ofthe Riccati equation (18).

ΔTHcmd(k)=−D(M)·G ^(T) ·P(M−1)ΦX′(k)+D(M)G ^(T) [C ^(T), ξ^(T)(M−2)C^(T),ξ^(T)(M−2)ξ^(T)(M−3)C ^(T), . . . . . . , ξ^(T)(M−2)ξ^(T)(M−3) . .. ξ^(T)(0)C ^(T)]Gcyl_cmd(k)

[0094] where

Gcyl_cmd(k)=[Gcyl_cmd(k+1), Gcyl_cmd(k+2), Gcyl_cmd(k+3), . . . ,Gycl_cmd(k+M)]  (17)

P(M−j)=Φ^(T) P(M−(j+1)Φ−Φ^(T) P(M−(j+1))GD(M)G ^(T) P(M−(j+1))Φ+C ^(T) C  (18)

[0095] where

D(M−j)=[H+G ^(T) P(M−(j+1))G} ⁻¹

ξ(M−(j+1))=Φ−GD(M−j)G ^(T) P(M−(j+1))Φj=M−1, M−2, . . . , 1

[0096] By defining initial conditions as shown in the equation (19), Pand D are recursively obtained.

P(0)=C^(T)C

D(1)=[H+G ^(T) P(0)G] ⁻¹  (19)

[0097] In case of M=1 (when a one-step-ahead desired value isavailable), the equation (16) is expressed as shown by the equation (20)and the equation (17) is expressed as shown by the equation (21).

J _(G)={Gcyl_cmd(k+1)−Gcyl(k+1)}+HΔTHcmd(k)²  (20)

[0098] $\quad\begin{matrix}{{{\Delta \quad {{THcmd}(k)}} = {{{{- {D(1)}}G^{T}{P(0)}\quad \Phi \quad {X^{\prime}(k)}} + \quad {{D(1)}G^{T}C^{T}{Gcyl\_ cmd}\left( {k + 1} \right)}}\quad = {{{- \left\lbrack {H + {G^{T}{P(0)}G}} \right\rbrack^{- 1}}G^{T}{P(0)}\Phi \quad {X^{\prime}(k)}} + {\left\lbrack {H + {G^{T}{P(0)}G}} \right\rbrack^{- 1}G^{T}C^{T}{Gcyl\_ cmd}\left( {k + 1} \right)}}}}\quad} & (21)\end{matrix}$

[0099] The feedback coefficients for X′(k) and Gcyl_cmd(k+1) in theequation (21) are calculated. $\begin{matrix}{{{P(0)} = {{C^{T}C} = {{\begin{bmatrix}1 \\0 \\0 \\0\end{bmatrix}\begin{bmatrix}1 & 0 & 0 & 0\end{bmatrix}} = \begin{bmatrix}1 & 0 & 0 & 0 \\0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 \\0 & 0 & 0 & 0\end{bmatrix}}}}{{G^{T}{P(0)}G} = {{{\begin{bmatrix}0 & 0 & 0 & 1\end{bmatrix}\quad\begin{bmatrix}1 & 0 & 0 & 0 \\0 & 0 & 0 & 0 \\0 & 0 & 0 & 0 \\0 & 0 & 0 & 0\end{bmatrix}}\begin{bmatrix}0 \\0 \\0 \\1\end{bmatrix}} = {{\left. 0\quad\Downarrow \left\lbrack {H + {G^{T}{P(0)}G}} \right\rbrack^{- 1} \right.G^{T}{P(0)}\Phi} = {{{\begin{bmatrix}0 & 0 & 0 & 0\end{bmatrix}\left\lbrack {H + {G^{T}{P(0)}G}} \right\rbrack}^{- 1}G^{T}C^{T}} = {{\left. 0\quad\Downarrow \Delta \right.\quad {{THcmd}(k)}} = {{{- \begin{bmatrix}0 & 0 & 0 & 0\end{bmatrix}}{X^{\prime}(k)}} + {0{Gcyl\_ cmd}\left( {k + 1} \right)}}}}}}}} & (22)\end{matrix}$

[0100] Thus, when elements in the first row of the G and G′ vectors arezero due to the dead time, performing such conventional GPC does notimplement a feasible intake air amount feedback controller 23.

[0101] 1.3 Structure of Intake Air Amount FB Controller

[0102]FIG. 3 shows a structure of the intake air amount feedbackcontroller 23 in accordance with one embodiment of the presentinvention. According to the present invention, a feasible intake airamount feedback control using a preview control is implemented byconstructing the intake air amount FB controller 23 as shown in FIG. 3.

[0103] The intake air amount FB, controller 23 comprises an adaptivedisturbance observer 31, a state predictor 32 and a control unit 33. Theadaptive disturbance observer 31 identifies an estimated value γ1 fordisturbance that is applied to the intake manifold 3. The statepredictor 32 calculates a predicted value Pre_Gcyl for the output of theintake manifold 3, which is a plant, by using the estimated disturbanceγ1. Through a preview control algorithm using the predicted valuePre_Gcyl, the control unit 33 calculates a desired throttle openingangle THcmd, which is a control input to the plant. The control inputTHcmd includes a value obtained by multiplying the estimated disturbanceγ1 by a predetermined gain. The output Gcyl of the plant can converge toa desired value by causing the predicted value Pre_Gcyl to converge to adesired value. The above problem that the elements of the first row ofthe G and G′ vectors become zero is prevented by introducing the statepredictor 32. The state predictor 32 will be described.

[0104] Since a value required to compensate for the dead time dth causedby the electronic control throttle valve is Gcyl(k+dth−1), the modelequation (13) for the dynamic behavior of intake air is shifted by(dth−1) steps to the future.

Gcyl(k+dth−1)=Aair1·Gcyl(k+dth−2)+Aair2·Gcyl(k+dth−3)+Bair1·THcmd(k−1)  (23)

[0105] The equation (23) includes future values Gcyl(k+dth−2) andGcyl(k+dth−3) which cannot be observed. Therefore, these future valuesare erased. Such erasure may be achieved by recursive calculation asfollows. The equation (24) represents a prediction equation for theintake air amount Gcyl.

Gcyl(k+1)=Aair1·Gcyl(k)+Aair2·Gcyl(k−1)+Bair·THcmd(k−dth+1)  (a)

[0106]

Gcyl(k+2)=Aair1·Gcyl(k+1)+Aair2·Gcyl(k)+Bair·THcmd(k−dth+2)

[0107]

[0108] substitute the above equation (a) into Gcyl(k+1)

Gcyl(k+2)=Aair1·(Aair1·Gcyl(k)+Aair2·Gcyl(k−1)+Bair·THcmd(k−dth+1))+Aair2·Gcyl(k)++Bair·THcmd(k−dth+2)

[0109]

Gcyl(k+dth−1)=αair1·Gcyl(k)+αair2·Gcyl(k−1)++βair1·THcmd(k−1)+βair2·THcmd(k−2)+,. . . , +β_(airdth−1)·THcmd(k−dth+1)  (24)

[0110] where $\begin{bmatrix}{\alpha \quad {air}\quad 1} & {\alpha \quad {air}\quad 2} \\* & *\end{bmatrix} = \begin{bmatrix}{A\quad {air}\quad 1} & {A\quad {air}\quad 2} \\1 & 0\end{bmatrix}^{{dth} - 1}$

[0111] Assuming ${\begin{bmatrix}{A\quad {air}\quad 1} & {A\quad {air}\quad 2} \\1 & 0\end{bmatrix} = A},{\begin{bmatrix}{Bair1} \\0\end{bmatrix} = B},$

[0112] [βair1 βair2 . . . β_(airdth−2) β_(airdth−1)] is elements of thefirst row of [B AB . . . A^(dth−3)B A^(dth−4)B]

[0113] Although the GPC is a control theory using the principle ofoptimality, the GPC does not have sufficient robustness against modelingerrors and prediction errors because the GPC is not designed to considersuch errors. According to one embodiment of the invention, the estimateddisturbance γ1 is included in the prediction equation (24) so as tocompensate for modeling errors and prediction errors. The statepredictor calculates the equation (25) to determine the predicted valuePre_Gcyl.

Pre_Gcyl(k)=αair1·Gcyl(k)+αair2·Gcyl(k−1)+βair1·THcmd(k−1)+βair2·THcmd(k−2)+,. . . , +β_(airdth−1)·THcmd(k−dth+1)+γ1(k)≅Gcyl(k+dth−1)  (25)

[0114] Calculation of the predicted value by the state predictor 32allows the dead time to be compensated for, thereby enhancing a quickresponse of the intake air amount control. Since the estimateddisturbance γ1 is included in the predicted value, a steady-state errorbetween the output Gcyl of the intake manifold (which is an object to becontrolled) and the predicted value Pre_Gcyl can be eliminated.

[0115] The estimated disturbance γ1 is identified by the adaptivedisturbance observer 31. The adaptive disturbance observer 31 calculatesthe equation (26) to determine the estimated disturbance γ1.$\begin{matrix}{{{\gamma \quad 1(k)} = {{\gamma \quad 1\left( {k - 1} \right)} + {{\frac{{Pdov}\left( {k - 1} \right)}{1 + {{Pdov}\left( {k - 1} \right)}} \cdot {e\_ dov}}(k)}}}{{{where}\quad {e\_ dov}(k)} = {{{Gcyl}(k)} - {{Gcyl\_ hat}(k)}}}{{{{Gcyl\_ hat}(k)} = {{\alpha \quad {air}\quad {1 \cdot {{Gcyl}\left( {k - {dth} + 1} \right)}}} + {\alpha \quad {air}\quad {2 \cdot {{Gcyl}\left( {k - {dth}} \right)}}} + {\beta \quad {air}\quad {1 \cdot {{THcmd}\left( {k - {dth}} \right)}}} + {\beta \quad {air}\quad {2 \cdot {{THcmd}\left( {k - {dth} - 1} \right)}}} +}},\ldots \quad,{{{+ \beta_{{airdth} - 1}} \cdot {{THcmd}\left( {k - {2{dth}} + 2} \right)}} + {\gamma \quad 1\left( {k - 1} \right)}}}{{{Pdov}(k)} = {\frac{1}{\lambda_{1}}\left( {1 - \frac{\lambda_{2}{{Pdov}\left( {k - 1} \right)}}{\lambda_{1} + {\lambda_{2}{{Pdov}\left( {k - 1} \right)}}}} \right){{Pdov}\left( {k - 1} \right)}}}} & (26)\end{matrix}$

[0116] As apparent from the equation (26), the adaptive disturbanceobserver 31 calculates a predicted value Gcyl_hat(k) for a current cyclein a similar way to the prediction equation (25). The adaptivedisturbance observer 31 calculates an error e_dov between the predictedvalue Gcyl_hat(k) and the actually detected value Gcyl(k). A recursiveidentification algorithm is used to calculate the estimated disturbanceγ1 so that the error e_dov converges to zero. By using the recursiveidentification algorithm, the estimated disturbance can be quickly andsteadily identified. Further, even when noise is included in the outputof the intake manifold 3 which is a controlled object, variations in theestimated disturbance caused by such noise can be suppressed by theeffect of the statistical process by the recursive identificationalgorithm.

[0117] λ1 and λ2 are weighting parameters. The recursive identificationalgorithm is a least square method in case of λ1=1 and λ2=1, a weightedleast square method in case of λ1<1 and λ2=1, a fixed gain method incase of λ1=1 and λ2=0, and a gradually-decreasing gain method in case ofλ=1 and λ2<1.

[0118] Next, the control unit 33 will be described. The equation (27) isobtained by shifting the prediction equation (25) by one step to thefuture and then converting it to include future values. Such conversionto include the future values can be achieved by reversing the conversionprocess from the equation (23) to the equation (24). $\begin{matrix}{{{{Pre\_ Gcyl}(k)} = {{\alpha \quad {air}\quad {1 \cdot {Gcyl}}(k)} + {\alpha \quad {air}\quad {2 \cdot {Gcyl}}\left( {k - 1} \right)} + {\beta \quad {air}\quad {1 \cdot {{THcmd}\left( {k - 1} \right)}}} + {\beta \quad {air}\quad {2 \cdot {{THcmd}\left( {k - 2} \right)}}} +}},\ldots \quad,{{{{+ \beta_{{airdth} - 1}} \cdot {{THcmd}\left( {k - {dth} + 1} \right)}} + {\gamma \quad 1(k)}}\quad \cong \left. {{Gcyl}\left( {k + {dth} - 1} \right)}\Downarrow {{Gcyl}\left( {k + {dth}} \right)} \right. \cong {{\alpha \quad {air}\quad {1 \cdot {Gcyl}}\left( {k + 1} \right)} + {\alpha \quad {air}\quad {2 \cdot {Gcyl}}(k)} + {\beta \quad {air}\quad {1 \cdot {THcmd}}(k)\beta \quad {air}\quad {2 \cdot {{THcmd}\left( {k - 1} \right)}}} +}},\ldots \quad,{{{\beta_{{airdth} - 1} \cdot {{THcmd}\left( {k - {dth} + 2} \right)}} + {\gamma \quad 1\left. \left( {k + 1} \right)\Downarrow {{Gcyl}\left( {k + {dth}} \right)} \right.}} \cong {{{Aair1} \cdot {{Gcyl}\left( {k + {dth} - 1} \right)}} + {{Aair2} \cdot {{Gcyl}\left( {k + {dth} - 2} \right)}} + {B\quad {{air1} \cdot {{THcmd}(k)}}} + {{\gamma 1}\left( {k + 1} \right)}}}} & (27)\end{matrix}$

[0119] A difference operator Δ that is defined as Δ=1−Z⁻¹ is introducedto define an augmented system as shown by the equation (28). Egc is anerror between the actual intake air amount Gcyl and a desired valueGcyl_cmd for the intake air amount. One additional row is specified inthe augmented system so as to derive an integral term that suppresses asteady-state error. It is assumed that a variation in the disturbance isconstant (that is, Δγ1(k)=Δγ1(k+1)). $\begin{matrix}{{{{Gcyl}\left( {k + {dth}} \right)} = {{{{Aair}\quad {1 \cdot {{Gcyl}\left( {k + {dth} - 1} \right)}}} + {{Aair}\quad {2 \cdot {{Gcyl}\left( {k + {dth} - 2} \right)}}} + {B\quad {air}\quad {1 \cdot {{THcmd}(k)}}} + {{{\gamma 1}\left( {k + 1} \right)}{{Egc}(k)}}} = {{{{Gcyl}(k)} - {{Gcyl\_ cmd}\left. (k)\quad \quad\Downarrow \Delta \right.\quad {{Gcyl}\left( {k + {dth}} \right)}}} = {{{Aair}\quad {1 \cdot \Delta}\quad {{Gcyl}\left( {k + {dth} - 1} \right)}} + {{Aair}\quad {2 \cdot \Delta}\quad {{Gcyl}\left( {k + {dth} - 2} \right)}} + \quad {{Bair}\quad {1 \cdot \Delta}\quad {{THcmd}(k)}} + {{\Delta\gamma 1}\left( {k + 1} \right)}}}}}\quad {{{Egc}(k)} = {{{{Ehc}\left( {k - 1} \right)} + {\Delta \quad {{Egc}(k)}}} = {{{{Egc}\left( {k - 1} \right)} + {\Delta \quad {{Gcyl}(k)}} - {\Delta \quad {Gcyl\_ cmd}(k)}} = \quad {{{{Egc}\left( {k - 1} \right)} + {{Aair}\quad {1 \cdot \Delta}\quad {{Gcyl}\left( {k - 1} \right)}} + {{Aair}\quad {2 \cdot \Delta}\quad {{Gcyl}\left( {k - 2} \right)}} + \quad \quad {{Bair}\quad {1 \cdot \Delta}\quad {{THcdm}\left( {k - {dth}} \right)}} + {\Delta \quad \gamma \quad 1\left( {k - {dth} + 1} \right)} - {\Delta \quad {Gcyl\_ cmd}\left. (k)\quad\Downarrow\quad \begin{bmatrix}{{Egc}\left( {k + {dth}} \right)} \\{\Delta \quad {{Gcyl}\left( {k + {dth}} \right)}} \\{\Delta \quad {{Gcyl}\left( {k + {dth} - 1} \right)}}\end{bmatrix} \right.}} = {\begin{bmatrix}1 & {A\quad {air}\quad 1} & {A\quad {air}\quad 2} \\0 & {A\quad {air}\quad 1} & {A\quad {air}\quad 2} \\0 & 1 & 0\end{bmatrix}{\quad{{\begin{bmatrix}{{Egc}\left( {k + {dth} - 1} \right)} \\{\Delta \quad {{Gcyl}\left( {k + {dth} - 1} \right)}} \\{\Delta \quad {{Gcyl}\left( {k + {dth} - 2} \right)}}\end{bmatrix} + {\begin{bmatrix}{B\quad {air}\quad 1} \\{B\quad {air}\quad 1} \\0\end{bmatrix}\Delta \quad {{THCmd}(k)}} + {\begin{bmatrix}1 \\1 \\0\end{bmatrix}\Delta \quad \gamma \quad 1\left( {k + 1} \right)} + {\begin{bmatrix}{- 1} \\0 \\0\end{bmatrix}\Delta \quad {Gcyl\_ cmd}\left. \left( {k + {dth}} \right)\quad\Downarrow\Delta \right.\quad \gamma \quad 1\left( {k + 1} \right)}} = {{{\Delta\gamma}\quad 1{(k)\begin{bmatrix}{{Egc}\left( {k + {dth}} \right)} \\{\Delta \quad {{Gcyl}\left( {k + {dth}} \right)}} \\{\Delta \quad {{Gcyl}\left( {k + {dth} - 1} \right)}}\end{bmatrix}}} = {\begin{bmatrix}1 & {A\quad {air}\quad 1} & {A\quad {air}\quad 2} \\0 & {A\quad {air}\quad 1} & {A\quad {air}\quad 2} \\0 & 1 & 0\end{bmatrix}{\quad{{\begin{bmatrix}{{Egc}\left( {k + {dth} - 1} \right)} \\{\Delta \quad {{Gcyl}\left( {k + {dth} - 1} \right)}} \\{\Delta \quad {{Gcyl}\left( {k + {dth} - 2} \right)}}\end{bmatrix} + {\begin{bmatrix}{B\quad {air}\quad 1} \\{B\quad {air}\quad 1} \\0\end{bmatrix}\Delta \quad {{THCmd}(k)}} + {\begin{bmatrix}1 \\1 \\0\end{bmatrix}\Delta \quad \gamma \quad 1(k)} + {\begin{bmatrix}{- 1} \\0 \\0\end{bmatrix}\Delta \quad {Gcyl\_ cmd}\left. \left( {k + {dth}} \right)\quad\Downarrow {X\left( {k + 1} \right)} \right.}} = {{{\Psi \quad {X(k)}} + {\Gamma \quad \Delta \quad {{ThCmd}(k)}} + {\Omega \quad \Delta \quad \gamma \quad 1(k)} + {{\Theta\Delta}\quad {Gcyl\_ cmd}\left( {k + {dth}} \right){X(k)}}} = {{\begin{bmatrix}{{Egc}\left( {k + {dth} - 1} \right)} \\{\Delta \quad {{Gcyl}\left( {k + {dth} - 1} \right)}} \\{\Delta \quad {{Gcyl}\left( {k + {dth} - 2} \right)}}\end{bmatrix}\quad \Psi} = {{\begin{bmatrix}1 & {A\quad {air}\quad 1} & {A\quad {air}\quad 2} \\0 & {A\quad {air}\quad 1} & {A\quad {air}\quad 2} \\0 & 1 & 0\end{bmatrix}\quad \Gamma} = {{\begin{bmatrix}{B\quad {air}\quad 1} \\{B\quad {air}\quad 1} \\0\end{bmatrix}\quad \Omega} = {{\begin{bmatrix}1 \\1 \\0\end{bmatrix}\quad \Theta} = \begin{bmatrix}{- 1} \\0 \\0\end{bmatrix}}}}}}}}}}}}}}}}}} & (28)\end{matrix}$

[0120] A cost function J_(SG) is defined. Assuming that the number ofdesired preview stages is represented by Nr and the number ofdisturbance preview stages is represented by Nd, a cost function JSG isdefined by using N=Max(Nr, Nd). The number of the desired preview stagesNr is equivalent to the above-described time period M, and specifies atime period during which future values for the desired value Gcyl_cmdare to be used. The number of the disturbance preview stages Ndspecifies a time period during which future values for the estimateddisturbance γ1, which is calculated by the adaptive disturbance observer31, are to be used. In the present embodiment, Nd is zero and Nr is one.Accordingly, N=1. The cost function J_(SG) is expressed as shown by theequation (29). $\begin{matrix}{{J_{SG} = {{\sum\limits_{j = 1}^{N}\quad \left. \left\lbrack {{{X^{T}\left( {k + j} \right)} \cdot Q \cdot {X\left( {k + j} \right)}} + {{R \cdot \Delta}\quad {{THcmd}\left( {k + j - 1} \right)}^{2}}} \right\rbrack \quad\Downarrow N \right.} = 1}}{J_{SG} = {{{X^{T}\left( {k + 1} \right)} \cdot Q \cdot {X\left( {k + 1} \right)}} + {{R \cdot \Delta}\quad {{THcmd}(k)}^{2}}}}} & (29)\end{matrix}$

[0121] where Q is weighting parameter of Quantity of State X$Q = {\begin{bmatrix}{q1} & 0 & 0 \\0 & {q2} & 0 \\0 & 0 & 0\end{bmatrix}\quad \left( {{q1},{q2},{> 0}} \right)}$

[0122] R is input weighting parameter

[0123] A control input ΔTHcmd that minimizes the cost function J_(SG) isobtained by using the principle of optimality. When the solution Π ofthe Riccati equation shown in the equation (31) is used, the controlinput ΔTHcmd is expressed as shown by the equation (30). $\begin{matrix}{{\Delta \quad {{THcmd}(k)}} = {{{Fx} \cdot {X(k)}} + {\sum\limits_{i = 1}^{Nr}\quad {{{{Fr}(j)} \cdot \Delta}\quad {Gcyl\_ cmd}\left( {k + {dth} - 1 + 1} \right)}} + {\sum\limits_{i = 0}^{Nd}\quad {{{{Fd}(1)} \cdot \Delta}\quad \gamma \quad 1\left( {k + 1} \right)}}}} & (30)\end{matrix}$

[0124] where: Fx = −Λ(N)Γ^(T)Π(N − 1)Ψ${{Fr}(i)} = \left\{ {{\begin{matrix}{{{- {\Lambda (N)}}\Gamma^{T}{\Pi \left( {N - 1} \right)}\Theta}\quad} & \left( {i = 1} \right) \\{{- {\Lambda (N)}}\Gamma^{T}\left\{ {\prod\limits_{m = 2}^{i}\quad {{\zeta^{T}\left( {N - 2} \right)}{\zeta^{T}\left( {N - 3} \right)}\quad \cdots \quad {\zeta^{T}\left( {N - i} \right)}}} \right\} {\Pi \left( {N - 1} \right)}\Theta} & \left( {1 \geq 2} \right)\end{matrix}{{Fd}(1)}} = \left\{ \begin{matrix}{{{- {\Lambda (N)}}\Gamma^{T}{\Pi \left( {N - 1} \right)}\Omega}\quad} & \left( {1 = 0} \right) \\{{- {\Lambda (N)}}\Gamma^{T}\left\{ {\prod\limits_{m = 1}^{1}{{\zeta^{T}\left( {N - 1} \right)}{\zeta^{T}\left( {N - 2} \right)}\quad {\zeta^{T}\left( {N - 3} \right)}\quad {{\cdots\zeta}^{T}\left( {N - 1} \right)}}} \right\} {\Pi \left( {N - 1} \right)}\Omega} & \left( {1 \geq 1} \right)\end{matrix} \right.} \right.$

Π(N−j)=Q+Ψ^(T)Π(N−(j+1))Ψ−Ψ^(T)Π(N−(j+1))ΓΛ(N)Γ^(T)Π(N−(j+1))ΨΛ(N−j)=[R+Γ^(T)Π(N−(j+1))Γ]⁻¹ ζ(N−(j+1))=Ψ−ΓΛ(N−j)Γ^(T)Π(N−(j+1))Ψ j=N−1, N−2, . .. ,1  (31)

[0125] The equation (30) is solved based on the initial conditions ofthe equation (32). The equation (33) is derived in case of N=1 (that is,Nr=1 and Nd=0).

Π*0)=Q

Λ(1)=[R+Γ^(T)Π(0)Γ]⁻¹  (32)

ΔTHcmd(k)=Fx·X(k)+Fr(1)ΔGcyl_cmd(k+dth)+Fd(0)Δγ1(k)  (33)

[0126] where:

Fx=−Λ(1)Γ^(T)Π(0)Ψ

Fr(1)=−Λ(1)Γ^(T)Π(0)Θ

Fd(0)=−Λ(1)Γ^(T)Π(0)Ω

[0127] Feedback coefficients Fx, Fd and a feedforward coefficient Fr inthe equation (33) are calculated as follows. $\begin{matrix}{{{\Lambda (1)} = \frac{1}{R + {1B\quad {air}\quad 1^{2}} + {{q2Bair}\quad 1^{2}}}}{{Fx} = \left\lbrack {\frac{{- {q1Bair}}\quad 1}{R + {{q1Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}{Aair}\quad 1{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}{Aair}\quad 2{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}} \right\rbrack}{{{Fr}(1)} = \frac{q\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}}{{{Fd}(0)} = \frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}}} & (34)\end{matrix}$

[0128] The control input ΔTHcmd in case of N=1 is calculated based onthe equations (33) and (34).

ΔTHcmd(k)=Fx·X(k)+Fr(1)ΔGcyl_cmd(k+dth)+Fd(0)Δγ1(k)

[0129] $\begin{matrix}{{{\Delta \quad {{THcmd}(k)}} = {{{{Fx} \cdot {X(k)}} + {{{Fr}(1)}\Delta \quad {Gcyl\_ cmd}\left( {k + {dth}} \right)} + {{{Fd}(0)}\Delta \quad \gamma \quad 1(k)}} = {{\frac{{- q}\quad 1{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}{{Egc}\left( {k + {dth} - 1} \right)}} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 1\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}\Delta \quad {{Gcyl}\left( {k + {dth} - 1} \right)}} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 2\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}\Delta \quad {{Gcyl}\left( {k + {dth} - 2} \right)}} + {\frac{{q1Bair}\quad 1}{R + {{q1Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}\Delta \quad {Gcyl\_ cmd}\left( {k + {dth}} \right)} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}\Delta \quad \gamma \quad 1(k)}}}}\quad} & (35)\end{matrix}$

[0130] The equation (35) is an equation for calculating the differentialΔTHcmd. The control input THcmd is calculated by integrating theequation (35). $\begin{matrix}{{{{THcmd}(k)} = {{\frac{{- q}\quad 1{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}{\sum\limits_{i = 0}^{k}\quad {{Egc}\left( {i + {dth} - 1} \right)}}} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 1\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}{Gcyl}\left( {k + {dth} - 1} \right)} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 2\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}\quad {{Gcyl}\left( {k + {dth} - 2} \right)}} + {\frac{q\quad 1{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}\quad {Gcyl\_ cmd}\left( {k + {dth}} \right)} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}\quad \gamma \quad 1(k)}\quad - {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 1\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}{{Gcyl}\left( {0 + {dth} - 1} \right)}} - {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 1\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}{{Gcyl}\left( {0 + {dth} - 2} \right)}} - {\frac{q\quad 1{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}\quad {Gcyl\_ cmd}\left( {0 + {dth}} \right)} - {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}\quad \gamma \quad 1(0)}}}\quad} & (36)\end{matrix}$

[0131] Assuming that initial values of Gcyl(1+dth−1) through Gcyl(0),Gcyl_cmd(0+dth) through Gcyl_cmd(0), γ1(0) and THcmd(0) are zero, theequation (36) is expressed as shown by the equation (37).$\begin{matrix}{{{{THcmd}(k)} = {{\frac{{- q}\quad 1{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}{\sum\limits_{i = 0}^{k}\quad {{Egc}\left( {i + {dth} - 1} \right)}}} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 1\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}{Gcyl}\left( {k + {dth} - 1} \right)} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 2\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}\quad {{Gcyl}\left( {k + {dth} - 2} \right)}} + {\frac{q\quad 1{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}\quad {Gcyl\_ cmd}\left( {k + {dth}} \right)} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}\quad \gamma \quad 1(k)}}}\quad} & (37)\end{matrix}$

[0132] The equation (37) contains future values Gcyl(k+dth−1) andGcyl(k+dth−2) that cannot be observed at the current time point “k.”Instead of these values, predicted values Pre_Gcyl(k) and Pre_Gcyl(k−1)calculated by the state predictor 32 are used. The equation (38) isexecuted by the control unit 33. Thus, the control input THcmd (k) isgenerated by the control unit 33. $\begin{matrix}{{{{{THcmd}(k)}=={{\frac{{- q}\quad 1{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}{\sum\limits_{i = 0}^{k}\quad {{Pre\_ Egc}(k)}}} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 1\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}{Pre\_ Gcyl}(k)} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 2\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}{Pre\_ Gcyl}\left( {k - 1} \right)} + {\frac{q\quad 1{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}\quad {Gcyl\_ cmd}\left( {k + {dth}} \right)} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}\quad \gamma \quad 1(k)\quad {where}\quad {Pre\_ Egc}(k)}}} = {{{{Pre\_ Gcyl}(k)} - {{Gcyl\_ cmd}\left( {k + {dth} - 1} \right)}} \cong {{{Gcyl}\left( {k + {dth} - 1} \right)} - {{Gcyl\_ cmd}\left( {k + {dth} - 1} \right)}}}}\quad} & (38)\end{matrix}$

[0133] Since a feedback term of the estimated disturbance value γ1 iscontained in the control input THcmd, an error between the intake airamount Gcyl and the desired value Gcyl_cmd, which may be caused by theapplication of disturbance, can quickly converge. Since the feedforwardterm Gcyl_cmd(k+dth) for the desired value is contained in the controlinput THcmd, the capability that the intake air amount Gcyl follows thedesired value Gcyl_cmd is improved.

[0134] 1.4 Result of Simulation of Intake Air Amount FB Control

[0135]FIG. 4 shows a model for a virtual controlled object that is usedin a simulation of the intake air amount FB control in accordance withone embodiment of the present invention. The virtual controlled objecthas a structure that is based on the model equation (13). A controlinput is a desired throttle opening angle THcmd(k−dth) that is delayedby a time “dth”. A control output is an intake air amount Gcyl(k). Theintake air amount one cycle ago Gcyl(k−1) and the intake air amount twocycles ago Gcyl(k−2) are fed back.

[0136] The simulation is structured to add three disturbances to thevirtual controlled object. An input disturbance d1, a state-quantitydisturbance d2 and an output disturbance d3 are shown in FIG. 4. Theinput disturbance d1 includes, for example, a variation in behavior ofthe throttle valve. The state-quantity disturbance d2 includes, forexample, a modeling error. The output disturbance d3 includes, forexample, noise of sensors.

[0137] Table 1 shows conditions of case G-1 through G-5 performed in thesimulation. TABLE 1 use of the F/B use of the term for the use of theF/F term for estimated estimated the desired disturbance γ 1 disturbancevalue in case disturbance in the control γ 1 in state the control numberd1, d2, d3 unit 33 predictor 32 unit 33 G-1 ◯ ◯ ◯ G-2 ◯ ◯ G-3 ◯ ◯ ◯ G-4◯ ◯ ◯ ◯ G-5 ◯ ◯ ◯

[0138] In the case G-1, no disturbance is added. In the state predictor32 and the control unit 33, the predicted value Pre_Gcyl and the desiredthrottle opening angle THcmd are calculated by using the estimateddisturbance value γ1. FIG. 5 shows a result of the simulation case G-1.Since there exists no disturbance, there is no error between thepredicted value Pre_Gcyl and the actual intake air amount Gcyl. Thecontrol unit 33 can cause the intake air amount Gcyl to follow thedesired value Gcyl_cmd without generating any overshooting.

[0139] In the case G-2, the disturbances d1 through d3 are added and theestimated disturbance value γ1 is not used in either the state predictor32 or the control unit 33. FIG. 6 shows a result of the simulation caseG-2. A steady-state error is generated between the predicted valuePre_Gcyl and the actual intake air amount Gcyl due to the disturbance.Since the control input THcmd is calculated based on the predicted valuePre_Gcyl, the control unit 33 cannot cause the intake air amount Gcyl toconverge to the desired value Gcyl_cmd.

[0140] In the case G-3, the disturbances d1 through d3 are added and theestimated disturbance value γ1 is used in the predictor 32. However, theestimated disturbance value γ1 is not used in the control unit 33. FIG.7 shows a result of the simulation case G-3. A steady-state errorbetween Pre_Gcyl and Gcyl, which is caused by the disturbances, isremoved by the predictor 32. Accordingly, the control unit 33 can causethe intake air amount Gcyl to converge to the desired value Gcyl_cmd.However, a convergence speed is relatively slow since a feedback termbased on the estimated disturbance value γ1 is not included in thecontrol input THcmd.

[0141] In the case G-4, the disturbances d1 through d3 are added and theestimated disturbance value γ1 is used in both of the predictor 32 andthe control unit 33. The case G-4 corresponds to a preferable embodimentaccording to the present invention that has been described abovereferring to FIG. 3. FIG. 8 shows a result of the simulation case G-4.As apparent from comparison with FIG. 7, the case G-4 significantlydecreases the time required for the error between the intake air amountGcyl and the desired value Gcyl_cmd to converge.

[0142] In the case G-5, the disturbances d1 through d3 are added and theestimated disturbance value γ1 is used in both of the predictor 32 andthe control unit 33. However, the feedforward term Gcyl_cmd(k+dth) ofthe desired value is not included in the control input THcmd. FIG. 9shows a result of the simulation case G-5. As apparent from comparisonwith FIG. 8, a speed that the intake air amount Gcyl follows the desiredvalue Gcyl_cmd slows down. This is because it is only an integral term(i.e., Pre_Egc term) that the control input THcmd includes for the errorbetween Gcyl and Gcyl_cmd. Thus, the capability that the intake airamount Gcyl follows the desired value Gcyl_cmd can be improved byincluding a feedforward term for the desired value in the control input.

[0143] Here, a case in which the control model has no dead time will bestudied. In such a case, the state predictor 32 may be removed. Themodel for the dynamic behavior of intake air for controlling the intakeair amount Gcyl can be expressed as shown by the equation (39).

Gcyl(k+1)=Aair1·Gcyl(k)+Aair2·Gcyl(k−1)+Bair1·THcmd(k)  (39)

[0144] Since there exists no dead time, the equation (26) performed bythe adaptive disturbance observer 31 is expressed by the equation (40).$\begin{matrix}{\quad {{{\gamma \quad 1(k)} = {{\gamma \quad 1\left( {k - 1} \right)} + {\frac{{Pdov}\left( {k - 1} \right)}{1 + {{Pdov}\left( {k - 1} \right)}}{e\_ dov}(k)}}}{{{e\_ dov}(k)} = {{{Gcyl}(k)} - {{Gcyl\_ hat}(k)}}}\quad {{{Gcyl\_ hat}(k)} = {{{Aair}\quad {1 \cdot {{Gcyl}\left( {k - 1} \right)}}} + {{Aair}\quad {2 \cdot {{Gcyl}\left( {k - 2} \right)}}} + \quad {{Bair}\quad {1 \cdot {{THcmd}\left( {k - 1} \right)}}} + {\gamma \quad 1\left( {k - 1} \right)}}}\quad {{{Pdov}(k)} = {\frac{1}{\lambda_{1}}\left( {1 - \frac{\lambda_{2}{{Pdov}\left( {k - 1} \right)}}{\lambda_{1} + {\lambda_{2}{{Pdov}\left( {k - 1} \right)}}}} \right){{Pdov}\left( {k - 1} \right)}}}}} & (40)\end{matrix}$

[0145] Since there exists no dead time, the equation (38) performed bythe control unit 33 is expressed by the equation (41). $\begin{matrix}{{{{THcmd}(k)}=={{\frac{{- q}\quad 1{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}{\sum\limits_{i = 0}^{k}\quad {{Egc}(k)}}} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 1\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}{{Gcyl}(k)}} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}A\quad {air}\quad 2\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}{{Gcyl}\left( {k - 1} \right)}} + {\frac{q\quad 1{Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {{q2}\quad {Bair}\quad 1^{2}}}{Gcyl\_ cmd}\left( {k + {dth}} \right)} + {\frac{{- \left( {{q\quad 1} + {q\quad 2}} \right)}\quad {Bair}\quad 1}{R + {q\quad 1{Bair}\quad 1^{2}} + {q\quad 2{Bair}\quad 1^{2}}}\quad \gamma \quad 1(k)\quad {{Egc}(k)}}}} = {{{Gcyl}(k)} - {{Gcyl\_ cmd}(k)}}} & (41)\end{matrix}$

[0146] As to such a case that includes no dead time, a simulation shownin Table 2 has been performed. TABLE 2 use of the F/B term use of theF/F term for the estimated for the desired disturbance disturbance γ 1in value in the control case number d1, d2, d3 the control unit 33 unit33 G-6 ◯ ◯ ◯ G-7 ◯ ◯

[0147] In the case G-6, the disturbance d1 through d3 are added and thecontrol unit 33 calculates the desired throttle opening angle THcmd byusing the estimated disturbance value γ1. The control input THcmdincludes the feedforward term Gcyl_cmd(k+dth) of the desired value. FIG.10 shows a result of the simulation case G-6.

[0148] In the case G-7, the disturbances d1 through d3 are added and thecontrol unit 33 calculates the desired throttle opening angle THcmdwithout using the estimated disturbance value γ1. The control inputTHcmd includes the feedforward term Gcyl_cmd(k+dth) of the desiredvalue. FIG. 11 shows a result of the simulation case G-7.

[0149]FIG. 12 shows comparison between the behavior of Gcyl shown inFIG. 10 (case G-6) and the behavior of Gcyl shown in FIG. 11 (case G-7).It is seen that the former is better than the latter from a viewpoint ofconvergence characteristics. Thus, the convergence characteristics ofthe controlled variable Gcyl relative to the desired value Gcyl_cmd canbe improved by calculating the control input THcmd through the use ofthe estimated disturbance value γ1 in the control unit 33.

[0150] In the embodiments as described above, the adaptive disturbanceobserver using the recursive identification algorithm is used so as toestimate a disturbance. Alternatively, another appropriate estimator,which may estimate a disturbance referring to a predetermined map or thelike, may be used. Further, in the embodiments as described above, thethrottle valve is used as a valve for controlling the intake air amount.Alternatively, another valve capable of controlling the intake airamount, for example, a bypass valve, may be used.

[0151] 2. Rotational Speed Feedback Control

[0152] 2.1 Modeling of Engine

[0153] A scheme of modeling the engine 2 will be described. The engine 2is represented by a model in which its input is the intake air amountGcyl and its output is the rotational speed NE.

[0154] The equation of motion in an inertial system of the engine isexpressed by the equation (42). Here, Ieng denotes inertial moment(kgm²) of the engine, Kne denotes a friction coefficient of the engine,and NE denotes an engine rotational speed (rad/sec). Teng denotes atorque (Nm) of the engine, Tload denotes an equipment-driving torque(Nm) for driving electrical components such as an air-conditioner, apower generator and the like which are mounted on the vehicle. Tdrvdenotes a vehicle-driving torque (Nm) for driving the vehicle, which isdistributed to a driving system of the vehicle. “t” indicates time.

Ieng·NE(t)=−Kne·NE(t)+Teng(t)−Tload(t)−Tdrv(t)  (42)

[0155] The engine torque Teng is expressed as shown by the equation(43). Ktrq denotes a torque coefficient, which is determined inaccordance with the engine rotational speed NE, an ignition timing IG ofthe engine and an equivalence ratio λ (a reciprocal of the air/fuelratio). Gcyl denotes the amount of air (g) that is introduced into theengine.

Teng(t)=Ktrq·Gcyl(t)  (43)

[0156] The equation (44) is obtained by substituting the equation (43)into the equation (42). The equation (44) represents a first-order delaysystem for the rotational speed NE in which its input is the intake airamount Gcyl. “−(Tload+Tdrv)/Ieng” is added as a disturbance term.$\begin{matrix}{{\overset{.}{NE}(t)} = {{{- \frac{Kne}{Ieng}} \cdot {{NE}(t)}} + {\frac{Ktrq}{Ieng} \cdot {{Gcyl}(t)}} + {\frac{- 1}{Ieng}\left( {{{Tload}(t)} + {{Tdrv}(t)}} \right)}}} & (44)\end{matrix}$

[0157] The equation (44) is converted into a discrete-time system toderive the equation (45). “T” denotes the length of the sampling cycle.Each cycle is identified by “k”. The equation (45) is a model equationfor the inertial system of the engine.

NE(k+1)=Ane·NE(k)+Bne·Gcyl(k)+Cne·(Tload(k)+Tdrv(k))  (45)

[0158] where ${Ane} = ^{{- \frac{Kne}{Ieng}}T}$${Bne} = {\int_{0}^{T}{\frac{Ktrq}{Ieng}^{{- \frac{Kne}{Ieng}}{({T - \tau})}}\quad {\tau}}}$${Cne} = {\int_{0}^{T}{\frac{- 1}{Ieng}^{{- \frac{Kne}{Ieng}}{({T - \tau})}}\quad {\tau}}}$

[0159] Model parameters Ane, Bne and Cne vary in accordance with therotational speed NE and the throttle opening angle TH. The modelparameters based on the rotational speed NE and the throttle openingangle TH may be pre-stored in the memory 1 c as a map. Alternatively, anidentifier may be provided in the control unit to identify these modelparameters.110

[0160] 2.2 Structure of Rotational Sipeed FB Control Unit

[0161]FIG. 13 shows a block diagram of the rotational speed FBcontroller 21 in accordance with one embodiment of the presentinvention. The rotational speed FB controller 21 comprises an adaptivedisturbance observer 41, a state predictor 42 and a control unit 43. Theadaptive disturbance observer 41 and the state predictor 42 have asimilar structure to those shown in FIG. 3 for the intake air amount FBcontroller 21.

[0162] The adaptive disturbance observer 41 identifies an estimatedvalue 6 ne for disturbance that is applied to the engine 2. The statepredictor 42 calculates a predicted value Pre_NE for the output (i.e.,an engine rotational speed) of the engine, which is a plant, based onthe estimated disturbance δne. Through a response assignment controlusing a predicted value Pre_NE, the control unit 43 calculates a desiredintake air amount Gcyl_cmd, which is a control input to the plant. Thecontrol input Gcyl_cmd includes a value obtained by multiplying theestimated disturbance 8 ne by a predetermined gain. The output NE of theplant can converge to a desired value by causing the predicted valuePre_NE to converge to a desired value.

[0163] The state predictor 42 will be described. As described abovereferring to FIG. 2, the rotational speed FB controller 21 is disposedupstream of the intake air amount FB controller 23. A dead time isincluded in the dynamic behavior of intake air of the intake manifold 3.If this dead time is compensated for by both of the rotational speed FBcontrol and the intake air amount FB control, some interference mayoccur. Therefore, the dead time included in the intake manifold 3 iscompensated for by the intake air amount FB controller 23. The intakemanifold 3 is regarded as a dead time element from the rotational speedFB controller 21. As a result, Gcyl_cmd(k−dth)=Gcyl(k) is seen from therotational speed FB controller 21. In other words, the rotational speedFB control unit 21 sees that the intake air amount Gcyl is introducedinto the engine when a dead time dth elapses after the calculation ofGcyl_cmd. Accordingly, the model equation (45) of the engine inertialsystem can be expressed as shown by the equation (46). Here, adisturbance Td indicates a sum of Tload and Tdrv.

NE(k+1)=Ane·Ne(k)+Bne·Gcyl_cmd(k−dth)+Cne·Td(k)  (46)

[0164] s

[0165] In order to compensate for the dead time dth, a control outputNE(k+dth) needs to be predicted. The equation (46) is shifted by (dth-1)steps to the future.

NE(k+dth)=Ane·Ne(k+dth−1)+Bne·Gcyl_cmd(k−1)+Cne·Td(k+dth−1)  (47)

[0166] Since the equation (47) includes future values NE(k+dth−1) andTd(k+dth−1) that cannot be observed, these future values are erased.Such erasure can be performed in a similar way to the erasure of thefuture values from the equation (23) as described above.

NE(k+1)=Ane·NE(k)+Bne·Gcyl_cmd(k−dth)+Cne·Td(k)  (a)

[0167] NE(k +2) =Ane NE(k +1) +Bne Gcyl_ cmd(k - dth +1) +Cne * Td(k +1)

NE(k+2)=Ane·Ane·NE(k)+Bne·Gcyl_cmd(k−dth)+Cne·Td(k)

[0168]

NE(k+1)

(a)

NE(k+2)=Ane{Ane·NE(k)+Bne·Gcyl_cmd(k−dth)+Cne·Td(k)}+Bne·Gcyl_cmd(k+dth+1)+Cne·Td(k+1)

[0169]

NE(k+dth)=Ane^(dth)NE(k)+Bne Gcy_cmd(k−1)+Ane·Bne·Gcy_cmd(k−2)+ . . .+Ane^(dth−1)Bne·Gcyl_cmd(k−dth)+Cne·Td(k−dth−1)+Ane Cne·Td(k+dth−2)+ . .. +Ane^(dth−1)Cne·Td(k)  (48)

[0170] It is hard to predict Td(k+dth−1) through Td(k) in the equation(48) since they change in accordance with operations of the driverand/or traveling conditions. Therefore, it is assumed that thedisturbance Td is constant as shown by the equation (49). According tothis assumption, the equation (48) is expressed by the equation (50).

Td(k+dth−1)=Td(k+dth−2)= . . . =Td(k+1)=Td(k)  (49)

NE(k+dth)=Ane^(dth)·NE(k)+Bne·Gcyl_cmd(k−1)+Ane1]Bne·Gcyl_cmd(k−2)+ . .. +Ane^(dth−1)·Bne·Gcyl_cmd(k−dth)+(1+Ane+ . . .+Ane^(dth−1))Cne·Td(k)  (50)

[0171] An estimated disturbance value δne is introduced into theequation (60). The estimated disturbance value δne includes not only anestimation error of the disturbance Td but also other disturbancesapplied to the plant. The equation (51) is executed by the stateestimator 42 to calculate a predicted value Pre_NE of the rotationalspeed.

NE(k+dth)≅αnel·NE(k)+βnel·Gcyl_cmd(k−1)+βne2·Gcyl_cmd(k−2)+ . . . . . .+βnedth·Gcyl_cmd(k−dth)+γyne·(Td(k)+δne(k))=Pre_NE(k)  (51)

[0172] where αne1=Ane^(dth)

[0173] βnej=Ane^(j−1)Bne j=1,2, . . . , dth

[0174] γne=(1+Ane+ . . . +Ane^(dth−1))Cne

[0175] By determining the predicted value by the state predictor 42, thedead time is compensated for and hence a quick response of therotational speed control can be enhanced. Since the predicted valuePre_NE is calculated based on the estimated disturbance δne, asteady-state error between the output NE of the engine (which is acontrolled object) and the predicted value Pre_NE is eliminated.

[0176] The estimated disturbance δne is identified by the adaptivedisturbance observer 41. The adaptive disturbance observer 41 executesthe equation (52) to determine the estimated disturbance δne.$\begin{matrix}{{\delta \quad {{ne}(k)}} = {{\delta \quad {{ne}\left( {k - 1} \right)}} + {{\frac{{Pdnde}\left( {k - 1} \right)}{1 + {{Pdne}\left( {k - 1} \right)}} \cdot {e\_ dne}}(k)}}} & (52)\end{matrix}$

[0177] where e_dne=NE(k)−NE_hat(k)

NE_hat(k)=αne1·NE(k−dth)+βne1·Gcyl_cmd(k−dth−1)+βne2·Gcyl_cmd(k−dth−2)+. . . . . . +βnedth·Gcyl_cmd(k−2dth)+γne·Td(k−dth)+γne·δne(k−1)

[0178]${{Pdne}(k)} = {\frac{1}{\lambda_{1}}\left( {1 - \frac{\lambda_{2}{{Pdne}\left( {k - 1} \right)}}{\lambda_{1} + {\lambda_{2}{{Pdne}\left( {k - 1} \right)}}}} \right){{Pdne}\left( {k - 1} \right)}}$

[0179] As apparent from the equation (52), the adaptive disturbanceobserver 41 calculates a predicted value NE-hat(k) for the current cycle(this is calculated by shifting the equation (51) by “dth” steps to thepast). Here, it is assumed that the estimated disturbance δne isconstant; that is, δne(k−dth)=δne(k−1). The adaptive disturbanceobserver 41 further calculates an error e_dne between the predictedvalue NE_hat(k) and the actually-detected value NE(k). Thereafter, arecursive identification algorithm is used to calculate the estimateddisturbance value δne so that the error e_dne converges to zero.

[0180] By using the recursive identification algorithm, the estimateddisturbance δne can be quickly and stably estimated. As described above,λ1 and λ2 are weighting parameters, which are determined in accordancewith the type of the recursive identification algorithm.

[0181] Next, the control unit 43 will be described. The equation (53) isobtained by shifting the prediction expression (51) by one step to thefuture and then converting it to include future values. Such conversionto include the future values can be achieved by reversing the conversionprocess from the equation (47) to the equation (48). Here, it is assumedthat variations in the future values of the disturbance Td and theestimated disturbance value δne are constant, that is, Td(k+dth)=Td(k)and δne(k+dth)=δne(k). $\begin{matrix}{{{Pre\_ NE}(k)} = {{{{\alpha \quad {ne}\quad {1 \cdot {{NE}(k)}}} + {\beta \quad {ne}\quad {1 \cdot {Gcyl\_ cmd}}\left( {k - 1} \right)} + {\beta \quad {{ne2} \cdot {Gcyl\_ cmd}}\left( {k - 2} \right)} + \cdots + {\beta \quad {{nedth} \cdot {Gcyl\_ cmd}}\left( {k - {dth}} \right)} + {\gamma \quad {{ne}\left( {{{Td}(k)} + {\delta \quad {{ne}(k)}}} \right)}}} \cong \left. {{NE}\left( {k + {dth}} \right)}\quad\Downarrow {{Ne}\left( {k + {dth} + 1} \right)} \right.} = {{{\alpha \quad {ne}\quad {1 \cdot {{NE}\left( {k + 1} \right)}}} + {\beta \quad {ne}\quad {1 \cdot {Gcyl\_ cmd}}(k)} + {\beta \quad {ne}\quad {2 \cdot \quad {Gcyl\_ cmd}}\left( {k - 1} \right)} + \cdots + {\beta \quad {{nedth} \cdot {Gcyl\_ cmd}}\left( {k - {dth} + 1} \right)} + {\gamma \quad \left. {{ne}\left( {{{Td}\left( {k + 1} \right)} + {\delta \quad {{ne}\left( {k + 1} \right)}}} \right)}\quad\Downarrow {{NE}\left( {k + {dth} + 1} \right)} \right.}} = {{{{Ane} \cdot {{NE}\left( {k + {dth}} \right)}} + {{{Bne} \cdot {Gcyl\_ cmd}}(k)} + {{Cne} \cdot \left. \left( {{{Td}\left( {k + {dth}} \right)} + {\delta \quad {{ne}\left( {k + {dth}} \right)}}} \right)\quad\Downarrow {{NE}\left( {k + {dth} + 1} \right)} \right.}} = {{{Ane} \cdot {{NE}\left( {k + {dth}} \right)}} + {{{Bne} \cdot {Gcyl\_ cmd}}(k)} + {{Cne} \cdot \left( {{{Td}(k)} + {\delta \quad {{ne}(k)}}} \right)}}}}}} & (53)\end{matrix}$

[0182] A switching function δne is defined to perform a responseassignment control. The switching function a ne allows convergencebehavior of the actual rotational speed NE to a desired value NE_cmd forthe rotational speed to be specified. E_ne denotes an error between theactual rotational speed NE and the desired value NE_cmd.

σne(k)=E_ne(k)+S_ne·E_ne(k−1)  (54)

[0183] s  where E_ne(k)=NE(k)−NE_cmd(k)

[0184] A control input is determined so that the switching function σnebecomes zero.

σne(k)=0

E_ne(k)=−S_ne·E_ne(k−1)  (55)

[0185] The equation (55) represents a first-order delay system having noinput. In other words, the control unit 43 operates to confine thecontrolled variable E_ne onto such a first-order delay system as shownby the equation (55).

[0186]FIG. 14 shows a phase plane with E_ne(k) being the vertical axisand E_ne(k−1) being the horizontal axis. A switching line 61 expressedby the equation (55) is shown in the phase plane. Assuming that a point62 is an initial value of a state variable (E_ne(k−1), E_ne(k)), thecontrol unit 43 places the state variable onto the switching line 61 andthen confines it on the switching line 61. Since the state variable isconfined in the first-order delay system having no input, the statevariable can automatically converge to the origin (that is, E_ne(k),E_ne(k−1)=0) of the phase plane with time. By confining the statevariable onto the switching line 61, the state quantity can converge tothe origin without being affected by disturbances.

[0187] A setting parameter S_ne of the equation (55) is established tosatisfy −1<S_ne<1. It is preferable that the setting parameter is set tosatisfy −1<S_ne<0. This is because the first-order delay system of theequation (55) may become a vibration-stable system when S_ne has apositive value.

[0188] The setting parameter S_ne is a parameter for specifying aconvergence speed of the error E_ne. Referring to FIG. 15, graphs 63, 64and 65 show a convergence speed in cases of S_ne=−1, −0.8 and −0.5,respectively. As the absolute value of the setting parameter S_nebecomes smaller, the convergence speed of the, error E_ne becomesfaster.

[0189] The control unit 43 calculates a control input Upas in accordancewith the equation (56). An equivalent control input Ueq is an input forconfining the state variable onto the switching line. A reaching lawinput Urch is an input for placing the state variable on the switchingline.

Gcyl_cmd(k)=Upas(k)=Ueq(k)+Urch(k)  (56)

[0190] A method for determining the equivalent control input Ueq will bedescribed. The equivalent control input Ueq has a function of holdingthe state variable at a given position in the phase plane. Therefore, itis required to satisfy the equation (57).

σne(k+dth+1)=σne(k+dth)  (57)

[0191] Based on the equation (54), the equation (57) is expressed by theequation (58). $\begin{matrix}{{{{E\_ ne}\left( {k + {dth} + 1} \right)} + {{{S\_ ne} \cdot {E\_ ne}}\left( {k + {dth}} \right)}} = {{{{E\_ ne}\left( {k + {dth}} \right)} + {{{S\_ ne} \cdot {E\_ ne}}\left. \left( {k + {dth} - 1} \right)\quad\Downarrow {{NE}\left( {k + {dth} + 1} \right)} \right.} - {{NE\_ cmd}\left( {k + {dth} + 1} \right)} + {{S\_ ne}\left\{ {{{NE}\left( {k + {dth}} \right)} - {{NE\_ cmd}\left( {k + {dth}} \right)}} \right\}}} = {{{NE}\left( {k + {dth}} \right)} - {{NE\_ cmd}\left( {k + {dth}} \right)} + {{S\_ ne}\left\{ {{{NE}\left( {k + {dth} - 1} \right)} - {{NE\_ cmd}\left( {k + {dth} - 1} \right)}} \right\}}}}} & (58)\end{matrix}$

[0192] The equation (59) can be obtained by substituting the equation(53) into the equation (58). $\begin{matrix}{{{{Ane} \cdot {{NE}\left( {k + {dth}} \right)}} + {{{Bne} \cdot {Gcyl\_ cmd}}(k)} + {{Cne} \cdot {{Td}(k)}} + {{{Cne} \cdot \delta}\quad {{ne}(k)}} - {{NE\_ cmd}\left( {k + {dth} + 1} \right)} + {{S\_ ne} \cdot {{NE}\left( {k + {dth}} \right)}} - {{{S\_ ne} \cdot {NE\_ cmd}}\left( {k + {dth}} \right)}} = {{{NE}\left( {k + {dth}} \right)} - {{NE\_ cmd}\left( {k + {dth}} \right)} + {{S\_ ne} \cdot {{NE}\left( {k + {dth} - 1} \right)}} - {{{S\_ ne} \cdot {NE\_ cmd}}\left( {k + {dth} - 1} \right)}}} & (59)\end{matrix}$

[0193] The control input Ueq(k) is calculated by the equation (60).

Ueq(k)=1/Bne{(1−Ane−S_ne)NE(k+dth)+S ₁₃ ne·NE(k+dth−1)+NE_cmd(k+dth+1)+(S_ne−1)·NE_cmd(k+dth)−S_ne·NE_cmd(k+dth−1)−Cne·Td(k)−Cne·δne(k)}  (60)

[0194] The equation (60) includes future values NE(k+dth) andNE(k+dth−1) which cannot be observed at the current time point “k”.Instead of these values, predicted values Pre_NE(k) and Pre_NE(k−1)calculated by the state predictor 42 are used. The control unit 43executes the equation (61) to determine the equivalent control inputUeq(k).

Ueq(k)=1/Bne{(1−Ane−S_ne)Pre_NE(k)+S_ne·Pre_NE(k−1)+NE_cmd(k+dth+1)+(S_ne−1)·NE_cmd(k+dth)−S_ne·NE_cmd(k+dth−1)−Cne·Td(k)−Cne·δne(k)}  (61)

[0195] Thus, the equivalent control input Ueq includes a disturbancefeedback term δne and a disturbance feedforward term Td. Accordingly,the error between the rotational speed NE and the desired value NE_cmd,which may be caused by application of disturbances to the engine 2 orthe controlled object, can quickly converge.

[0196] The control unit 43 further executes the equation (62) todetermine a reaching law input Urch. F denotes a reaching law gain.$\begin{matrix}{{{Urch}(k)} = {\frac{- F}{Bne}\sigma \quad {{ne}\left( {k + {dth}} \right)}}} & (62)\end{matrix}$

[0197] where 0<F<2

[0198] 2.3 Result of Simulation of Rotational Speed FB Control

[0199]FIG. 16 shows a model of a virtual controlled object that is usedin a simulation of the rotational speed FB control in accordance withone embodiment of the present invention. The virtual controlled objecthas a structure that is based on the model equation (46) of the engine.A control input is a desired intake air amount Gcyl_cmd that is delayedby a time “dth”. A control output is a rotational speed NE. A drivingtorque Td is applied to the controlled object as a disturbance. Therotational speed NE one cycle ago is fed back.

[0200] The simulation is structured to add three disturbances to thevirtual controlled object. Three positions at which an input disturbanceL1, a state-quantity disturbance L2 and an output disturbance L3 are tobe applied are shown. The input disturbance L1 includes, for example, anestimation error for the driving torque Td. The state-quantitydisturbance L2 includes, for example, a modeling error. The outputdisturbance L3 includes, for example, noise of sensors. Table 3 showsconditions for cases N-1 through N-5 performed in the simulation. TABLE3 use of the F/B term for the use of the estimated F/F term for use ofthe Disturbance disturbance the driving adaptive law case number L1, L2,L3 δ ne torque Td input Uadp N-1 ◯ ◯ ◯ N-2 ◯ N-3 ◯ ◯ N-4 ◯ ◯ N-5 ◯ ◯ ◯ ◯

[0201] In the case N-1, the disturbances L1 through L3 are added. Theestimated disturbance value δne and the driving torque Td are used inboth of the predictor 42 and the control unit 43. The case N-1 is apreferable case based on the rotational speed FB control in FIG. 13 inaccordance with one embodiment of the present invention. FIG. 17 shows aresult of the simulation case N-1. Under the condition in which thedisturbances are applied, the rotational speed NE can converge to thedesired value NE_cmd without any steady-state error. A capability thatthe rotational speed NE follows the desired value NE_cmd when thedesired value NE_cmd changes is good.

[0202] In the case N-2, the estimated disturbance value δne and thedriving torque Td are not used in either the predictor 42 or the controlunit 43. FIG. 18 shows a result of the simulation case N-2. Asteady-state error occurs between the actual rotational speed NE and thepredicted value Pre_NE due to the disturbances. Since the control unit43 performs the response assignment control based on the predicted valuePre_NE , the control unit 43 cannot cause the actual rotational speed NEto converge to the desired value NE_cmd.

[0203] In the case N-3, the predictor 42 and the control unit 43 use thedriving torque Td. The estimated disturbance value δne is not used ineither the predictor 42 or the control unit 43. FIG. 19 shows a resultof the simulation case N-3. At time t1, an estimation error Td_error forthe driving torque Td is applied as a step input. In response, theestimated disturbance value δne increases. Since the feedforward termfor the driving torque Td is used, the error between the actualrotational speed NE and the predicted value Pre_NE keeps unchanged.

[0204] At time t2, the driving torque Td changes. The other disturbancesL2 and L3 are still applied to the plant. The error between the actualrotational speed NE and the predicted value Pre_NE does not change byvirtue of the feedforward term for the driving torque Td. However, sincethe estimated disturbance value δne is not used, a steady-state errorbetween the actual rotational speed NE and the predicted value Pre_NEcannot be eliminated, and hence the rotational speed NE cannot convergeto the desired value NE_cmd.

[0205] In the case N-4, the estimated disturbance value δne and thedriving torque Td are not used in either the predictor 42 or the controlunit 43. An adaptive law input Uadp is added to the control input in thecontrol unit 43. The adaptive law input Uadp is expressed by theequation (63). “G” denotes a gain of the adaptive law input.$\begin{matrix}{{{Uadp}(k)} = {\frac{- G}{Bne}{\sum\limits_{i - 0}^{k + {dth}}{\delta \quad {{ne}(i)}}}}} & (63)\end{matrix}$

[0206]FIG. 20 shows a result of the simulation case N-4. The predictedvalue Pre_NE converges to the desired value NE_cmd by the control unit43. However, since the predicted value Pre_NE is not calculated based onthe estimated disturbance value δne, a steady-state error between theactual rotational speed NE and the predicted value Pre_NE cannot beeliminated, and hence the rotational speed NE cannot converge to thedesired value NE_cmd.

[0207] In the case N-5, the driving torque Td and the estimateddisturbance value (δne are used in the predictor 42. In the control unit43, the driving torque Td is used. The control unit 43 does not use theestimated disturbance value δne. Instead, the adaptive law input Uadp isadded to the control input in the control unit 43.

[0208]FIG. 21 shows a result of the simulation case N-5. The convergencetime of the rotational speed NE can be shorten by increasing the gain Gof the adaptive law input Uadp. However, as apparent from comparisonwith FIG. 17, integral overshooting occurs when the gain G increases.

[0209] Now, a case having no dead time will be examined. In such a case,the state predictor may be removed. A model for an inertial system ofthe engine for controlling the engine rotational speed NE can beexpressed as shown by the equation (64).

NE(k+1)=Ane·NE(k)+Bne·Gcyl_cmd(k)+Cne·Td(k)  (64)

[0210] Since there exists no dead time, the equation (52), which isexecuted by the adaptive disturbance observer 41, is expressed by theequation (65). $\begin{matrix}{{\delta \quad {{ne}(k)}} = {{\delta \quad {{ne}\left( {k - 1} \right)}} + {{\frac{{Pdne}\left( {k - 1} \right)}{1 + {{Pdne}\left( {k - 1} \right)}} \cdot {e\_ dne}}(k)}}} & (65)\end{matrix}$

[0211] where e_dne=NE(k)−NE_hat(k)

NE_hat(k)=Ane1·NE(k−1)+Bne1·Gcyl_cmd(k−1)+δne·Td(k−1)+δne·δne(k−1)

[0212]${{Pdne}(k)} = {\frac{1}{\lambda_{1}}\left( {1 - \frac{\lambda_{2}{{Pdne}\left( {k - 1} \right)}}{\lambda_{1} + {\lambda_{2}{{Pdne}\left( {k - 1} \right)}}}} \right){{Pdne}\left( {k - 1} \right)}}$

[0213] Since there exists no dead time, the equations (61) and (62),which are executed by the control unit 43, are expressed by theequations (66) and (67), respectively.

Ueq(k)=1/Bne{(1−Ane−S_ne)NE(k)+S_ne·NE(k−1)+NE_cmd(k+1)+(S_ne−1)·NE_cmd(k)−S_ne−NE_cmd(k−1)−Cne·Td(k)−Cne·δne(k)}  (66)

[0214] $\begin{matrix}{{{Urch}(k)} = {\frac{- F}{Bne}\sigma \quad {{ne}(k)}}} & (67)\end{matrix}$

[0215] where 0<F<2

[0216] As to cases N-6 through N-9 having no dead time, simulations asshown in Table 4 have been carried out. TABLE 4 use of the F/B term forthe use of the F/F estimated term for the use of the case disturbancedisturbance driving adaptive law number L1, L2, L3 value δ ne torque Tdinput Uadp N-6 ◯ ◯ ◯ N-7 ◯ ◯ N-8 ◯ ◯ ◯ (G: smaller) N-9 ◯ ◯ ◯ (G:larger)

[0217] In the case N-6, the disturbances L1 through L3 are added. Theestimated disturbance value δne and the driving torque Td are used inthe control unit 43. The case N-6 is a preferable case based on therotational speed FB control in accordance with one embodiment of thepresent invention as described above. FIG. 22 shows a result of thesimulation case N-6. Under the condition in which the disturbances areapplied, the rotational speed NE can converge to the desired valueNE_cmd without generating any steady-state error. The capability thatthe rotational speed NE follows the desired value NE_cmd when thedesired value NE_cmd changes is good.

[0218] In the case N-7, the driving torque Td is used in the controlunit 43. The control unit 43 does not use the estimated disturbancevalue δne. FIG. 23 shows a result of the simulation case N-7. It is seenthat the error between the actual rotational speed NE and the desiredvalue NE_cmd increases each time the disturbance are applied. The actualrotational speed NE cannot converge to the desired value NE_cmd sincethe estimated disturbance value δne is not used.

[0219] In the case N-8, the estimated disturbance value δne is not usedin the control unit 43. The adaptive law input Uadp is added to thecontrol input. In the case N-8, a gain G having a relatively small valueis used. FIG. 24 shows a result of the simulation case N-8. In the caseN-9, a gain G having a relatively large value is used in the controlunit 43. FIG. 25 shows a result of the simulation case N-9. FIG. 26shows a comparison between the behavior of the rotational speed NE shownin FIG. 25 (case N-9) and the behavior of the rotational speed NE shownin FIG. 24 (case N-8). By increasing the gain G of the adaptive lawinput Uadp, the convergence time of the rotational speed NE is shorten.However, integral overshooting occurs when the gain G increases. Incontrast, by using the estimated disturbance value δne corresponding tothe case N-6, the rotational speed NE can converge to the desired valueNE_cmd without causing integral overshooting.

[0220] In the above-described embodiments, the adaptive disturbanceobserver using the recursive identification algorithm is used so as toestimate a disturbance. Alternatively, another appropriate estimatorthat refers to a predetermined map or the like may be used to determinea disturbance.

[0221] 3. Operation Flow

[0222]FIG. 27 shows a flowchart of the rotational speed FB control andthe intake air amount FB control in accordance with one embodiment ofthe present invention shown in FIG. 2. This flowchart can be applied toa vehicle using any of a manual transmission (MT), an automatic manualtransmission (automatic MT) and an automatic transmission (AT).

[0223] In step S1, it is determined whether the vehicle is idling orgear change is being carried out in the transmission. When the answer ofstep S1 is YES, the process proceeds to step S2 to perform theabove-described rotational speed FB control.

[0224] In step S2, a desired engine rotational speed NE_cmd isdetermined. For example, when the engine is idling, the desired valueNE_cmd is set to a value in accordance with traveling conditions,warming-up state and so on. When gear change is being carried out, thedesired value NE_cmd is set to a value in accordance with a vehiclespeed and a selected gear ratio.

[0225] In step S3, a map stored in the memory 1 c of the ECU 1 isreferred to based on the detected actual rotational speed NE to extractthe model parameters Ane, Bne and Cne. In step S4, the equipment-drivingtorque Tload and the vehicle-driving torque Tdrv are determined. Theequipment-driving torque Tload may be calculated in accordance with, forexample, an ON/OFF state of electrical components mounted on thevehicle. The vehicle-driving torque Tdrv may be calculated in accordancewith traveling resistance, clutch conditions and so on. Tload and Tdrvare added to determine a driving torque Td as a disturbance. In step S5,the above-described rotational speed FB control is carried out tocalculate the desired intake air amount Gcyl_cmd.

[0226] On the other hand, when the vehicle is in a normal runningcondition, the process proceeds to step S6, in which a desired enginetorque is determined. The desired engine torque may be calculated inaccordance with an opening angle of the accelerator pedal, vehiclespeed, selected gear ratio, running environment and so on. In step S7,the intake air amount Gcyl_cmd required for implementing the desiredengine torque is calculated. For example, the desired intake air amountGcyl_cmd can be determined by referring to a predetermined map based onthe air/fuel ratio and the ignition timing.

[0227] In step S8, the intake air amount Gcyl is estimated. The intakeair amount Gcyl can be estimated based on the outputs from the airflowmeter 10 and the Pb sensor 11. In step S9, a predetermined map isreferred to based on the engine rotational speed NE and the throttleopening angle TH to extract the model parameters Aair1, Aair2 and Bair1.Instead of an opening angle of the throttle valve, the air amount Gththat passes through the throttle valve or the intake air amount Gcylcalculated in step S8 may be used.

[0228] In step S10, the above-described intake air amount FB control iscarried out to calculate the desired throttle opening angle THcmd.

[0229] It should be noted that a control scheme according to theinvention may be applied to various objects. A preview control accordingto the invention may be applied to various objects. A responseassignment control according to the invention may be also applied tovarious objects.

[0230] The invention may be applied to an engine to be used in avessel-propelling machine such as an outboard motor in which acrankshaft is disposed in the perpendicular direction.

What is claimed is:
 1. A control system for controlling a modeled plant,comprising a controller configured to; estimate disturbance applied tothe plant; and determine an input to the plant so that an output of theplant converges to a desired value, the input determined to include avalue obtained by multiplying the estimated disturbance by apredetermined gain.
 2. The control system of claim 1, wherein thecontroller is further configured to use a preview control algorithm todetermine the input to the plant.
 3. The control system of claim 1,wherein the controller is further configured to use a responseassignment control algorithm to determine the input to the plant.
 4. Thecontrol system of claim 1, wherein the controller includes an adaptivedisturbance observer that uses a recursive identification algorithm toidentify the estimated disturbance.
 5. The control system of claim 2,wherein the controller is further configured to determine the input tothe plant to include a value obtained by multiplying by a predeterminedgain a desired value for the output of the plant.
 6. The control systemof claim 2, wherein the plant is an intake manifold connected to anengine, wherein the intake manifold is modeled so that an input of theplant is a desired value for an opening angle of a valve that controlsan amount of intake air into the intake manifold and an output of theplant is an amount of intake air into the engine.
 7. The control systemof claim 6, wherein the valve that controls an amount of intake air intothe intake manifold is a throttle valve provided in the intake manifold.8. The control system of claim 6, wherein the controller includes astorage device for storing model parameters for the modeled plant,wherein the controller is further configured to extract a modelparameter based on a detected engine rotational speed and a detectedopening angle of a throttle valve and to determine the input to theplant based on the extracted model parameter.
 9. The control system ofclaim 6, wherein the desired value for the output of the plant is adesired intake air amount, wherein the controller is further configuredto set in the desired intake air amount an amount of intake air requiredto implement a desired engine rotational speed when the engine is idlingor when gear change is being carried out, wherein the controller isfurther configured to set in the desired intake air amount an amount ofintake air required to implement a desired engine torque when the engineis in a normal running condition.
 10. The control system of claim 3,wherein the plant is an engine, wherein the engine is modeled so that aninput of the plant is a desired value for an amount of intake air intothe engine and an output of the plant is a rotational speed of theengine.
 11. The control system of claim 10, wherein the controllerincludes a storage device for storing model parameters for the modeledplant, wherein the controller is further configured to extract a modelparameter based on a detected engine rotational speed and to determinethe input to the plant based on the extracted model parameter.
 12. Thecontrol system of claim 10, wherein the controller is further configuredto determine the input to the plant to include a value obtained bymultiplying by a predetermined gain an estimated value for a torque thatis required to drive a vehicle on which the engine is mounted.
 13. Thecontrol system of claim 10, wherein the controller is further configuredto determine the input to the plant to include a value obtained bymultiplying by a predetermined gain an estimated value for a torque thatis required to drive an equipment on a vehicle on which the engine ismounted.
 14. The control system of claim 10, wherein the controller isfurther configured to perform the response assignment control algorithmwhen the engine is idling or when gear change is being carried out. 15.A control system for controlling a modeled plant, comprising acontroller configured to: estimate disturbance applied to the plant;predict an output of the plant based on the estimated disturbance anddead time included in the plant; and determine an input to the plant sothat the predicted output converges to a desired value for the output ofthe plant.
 16. The control system of claim 15, wherein the controller isfurther configured to use a preview control algorithm to determine theinput to the plant.
 17. The control system of claim 15, wherein thecontroller is further configured to use a response assignment controlalgorithm to determine the input to the plant.
 18. The control system ofclaim 15, wherein the controller is further configured to determine theinput to the plant to include a value obtained by multiplying theestimated disturbance by a predetermined gain.
 19. The control system ofclaim 15, wherein the controller includes an adaptive disturbanceobserver that uses a recursive identification algorithm to identify theestimated disturbance.
 20. The control system of claim 16, wherein thecontroller is further configured to determine the input to the plant toinclude a value obtained by multiplying a desired value for the outputof the plant by a predetermined gain.
 21. The control system of claim16, wherein the plant is an intake manifold connected to an engine,wherein the intake manifold is modeled so that an input of the plant isa desired value for an opening angle of a valve that controls an amountof intake air into the intake manifold and an output of the plant is anamount of intake air into the engine.
 22. The control system of claim21, wherein the controller includes a storage device for storing modelparameters for the modeled plant, wherein the controller is furtherconfigured to extract a model parameter based on a detected enginerotational speed and a detected opening angle of a throttle valve and todetermine the input to the plant based on the extracted model parameter.23. The control system of claim 17, wherein the plant is an engine,wherein the engine is modeled so that an input of the plant is a desiredvalue for an amount of intake air into the engine and an output of theplant is a rotational speed of the engine.
 24. The control system ofclaim 23, wherein the controller includes a storage device for storingmodel parameters for the modeled plant, wherein the controller isfurther configured to extract a model parameter based on a detectedengine rotational speed and to determine the input to the plant based onthe extracted model parameter.
 25. The control system of claim 23,wherein the controller is further configured to determine the input tothe plant to include a value obtained by multiplying by a predeterminedgain an estimated value for a torque that is required to drive a vehicleon which the engine is mounted.
 26. The control system of claim 23,wherein the controller is further configured to determine the input tothe plant to include a value obtained by multiplying by a predeterminedgain an estimated value for a torque that is required to drive anequipment on a vehicle on which the engine is mounted.
 27. The controlsystem of claim of claim 23, wherein the controller is furtherconfigured to predict the output of the plant based on an estimatedvalue for a torque that is required to drive a vehicle on which theengine is mounted.
 28. The control system of claim of claim 23, whereinthe controller is further configured to predict the output of the plantbased on an estimated value for a torque that is required to drive anequipment on a vehicle on which the engine is mounted.
 29. A method forcontrolling a modeled plant, comprising the steps of: (a) estimatingdisturbance applied to the plant; and (b) determining an input to theplant so that an output of the plant converges to a desired value, theinput determined to include a value obtained by multiplying theestimated disturbance by a predetermined gain.
 30. The method of claim29, wherein the step (b) further comprises the step of using a previewcontrol algorithm to determine the input to the plant.
 31. The method ofclaim 29, wherein the step (b) further comprises the step of using aresponse assignment control algorithm to determine the input to theplant.
 32. The method of claim 29, wherein the step (a) furthercomprises the step of using a recursive identification algorithm toidentify the estimated disturbance.
 33. The method of claim 30, whereinthe step (b) further comprises the step of determining the input to theplant to include a value obtained by multiplying a desired value for theoutput of the plant by a predetermined gain.
 34. The method of claim 30,wherein the plant is an intake manifold connected to an engine; whereinthe intake manifold is modeled so that an input of the plant is adesired value for an opening angle of a valve that controls an amount ofintake air into the intake manifold and an output of the plant is anamount of intake air into the engine.
 35. The method of claim 34,wherein the step (b) further comprises the steps of: determining a modelparameter based on a detected engine rotational speed and a detectedopening angle of a throttle valve; and determining the input to theplant based on the model parameter.
 36. The method of claim 34, furthercomprising the steps of: setting in the desired value for the output ofthe plant an amount of intake air required to implement a desired enginerotational speed when the engine is idling or when gear change is beingcarried out; and setting in the desired value for the output of theplant an amount of intake air required to implement a desired enginetorque when the engine is in a normal running condition.
 37. The methodof claim 31, wherein the plant is an engine; wherein the engine ismodeled so that an input of the plant is a desired value for an amountof intake air into the engine and an output of the plant is a rotationalspeed of the engine.
 38. The method of claim 37, wherein the step (b)further comprises the steps of: determining a model parameter based on adetected engine rotational speed; and determining the input to the plantbased on the model parameter.
 39. The method of claim 37, wherein thestep (b) further comprises the step determining the input to the plantto include a value obtained by multiplying by a predetermined gain anestimated value for a torque that is required to drive a vehicle onwhich the engine is mounted.
 40. The method of claim 37, wherein thestep (b) further comprises the step of: determining the input to theplant to include a value obtained by multiplying by a predetermined gainan estimated value for a torque that is required to drive an equipmenton a vehicle on which the engine is mounted.
 41. The method of claim 37,wherein the response assignment control algorithm is performed when theengine is idling or when gear change is being carried out.
 42. A methodfor controlling a modeled plant, comprising the steps of: (a) estimatingdisturbance applied to the plant; (b) predicting an output of the plantbased on the estimated disturbance and dead time included in the plant;and (c) determining an input to the plant so that the predicted outputconverges to a desired value for the output of the plant.
 43. The methodof claim 42, wherein the step (c) further comprises the step of using apreview control algorithm to determine the input to the plant.
 44. Themethod of claim 42, wherein the step (c) further comprises the step ofusing a response assignment control algorithm to determine the input tothe plant.
 45. The method of claim 42, wherein the step (c) furthercomprises the step of determining the input to the plant to include avalue obtained by multiplying the estimated disturbance by apredetermined gain.
 46. The method of claim 42, wherein the step (a)further comprises the step of using a recursive identification algorithmto identify the estimated disturbance.
 47. The method of claim 43,wherein the step (c) further comprises the step of determining the inputto the plant to include a vallue obtained by multiplying a desired valuefor the output of the plant by a predetermined gain.
 48. A controllerfor controlling a modeled plant, comprising: (a) means for estimatingdisturbance applied to the plant; and (b) means for determining an inputto the plant so that an output of the plant converges to a desiredvalue, the input determined to include a value obtained by multiplyingthe estimated disturbance by a predetermined gain.
 49. The controller ofclaim 48, wherein the means (b) further comprises means for using apreview control algorithm to determine the input to the plant.
 50. Thecontroller of claim 48, wherein the means (b) further comprises meansfor using a response assignment control algorithm to determine the inputto the plant.
 51. The controller of claim 48, wherein the means (a)further comprises means for using a recursive identification algorithmto identify the estimated disturbance.
 52. The controller of claim 49,wherein the means (b) further comprises means for determining the inputto the plant to include a value obtained by multiplying a desired valuefor the output of the plant by a predetermined gain.
 53. The controllerof claim 49, wherein the plant is an intake manifold connected to anengine, wherein the intake manifold is modeled so that an input of theplant is a desired value for an opening angle of a valve that controlsan amount of intake air into the intake manifold and an output of theplant is an amount of intake air into the engine.
 54. The controller ofclaim 53, wherein the means (b) further comprises: means for determininga model parameter based on a detected engine rotational speed and adetected opening angle of a throttle valve; and means for determiningthe input to the plant based on the model parameter.
 55. The controllerof claim 53, further comprising: means for setting in the desired valuefor the output of the plant an amount of intake air required toimplement a desired engine rotational speed when the engine is idling orwhen gear change is being carried out; and means for setting in thedesired value for the output of the plant an amount of intake airrequired to implement a desired engine torque when the engine is in anormal running condition.
 56. The controller of claim 50, wherein theplant is an engine, wherein the engine is modeled so that an input ofthe plant is a desired value for an amount of intake air into the engineand an output of the plant is a rotational speed of the engine.
 57. Thecontroller of claim 56, wherein the means (b) further comprises: meansfor determining a model parameter based on a detected engine rotationalspeed; and means for determining the input to the plant based on themodel parameter.
 58. The controller of claim 56, wherein the means (b)further comprises: means for determining the input to the plant toinclude a value obtained by multiplying by a predetermined gain anestimated value for a torque that is required to drive a vehicle onwhich the engine is mounted.
 59. The controller of claim 56, wherein themeans (b) further comprises: means for determining the input to theplant to include a value obtained by multiplying by a predetermined gainan estimated value for a torque that is required to drive an equipmenton a vehicle on which the engine is mounted.
 60. The controller of claim56, further comprising means for performing the response assignmentcontrol algorithm when the engine is idling or when gear change is beingcarried out.
 61. A controller for controlling a modeled plant,comprising: (a) means for estimating disturbance applied to the plant;(b) means for predicting an output of the plant based on the estimateddisturbance and dead time included in the plant; and (c) means fordetermining an input to the plant so that the predicted output convergesto a desired value for the output of the plant.
 62. The controller ofclaim 61, wherein the means (c) further comprises means for using apreview control algorithm to determine the input to the plant.
 63. Thecontroller of claim 61, wherein the means (c) further comprises meansfor using a response assignment control algorithm to determine the inputto the plant.
 64. The controller of claim 61, wherein the means (c)further comprises means for determining the input to the plant toinclude a value obtained by multiplying the estimated disturbance by apredetermined gain.
 65. The controller of claim 61, wherein the means(a) further comprises means for using a recursive identificationalgorithm to identify the estimated disturbance.
 66. The controller ofclaim 62, wherein the means (c) further comprises means for determiningthe input to the plant to include a value obtained by multiplying adesired value for the output of the plant by a predetermined gain.